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Portable  Neutron  Source 

A.  Modeling  and  simulation  tools: 

1.  Laser-target  deposition  model. 

A  two-dimensional  fully  relativistic  electromagnetic  particle-in-cell  model  (Ref.  [1])  was 
developed.  The  particle-in-cell  model  describes  the  evolution  of  the  target  by  moving  “quasi¬ 
particles”  representing  each  species  (electrons  and  ions).  The  “quasi-particles”  are  driven  by  the 
laser  electromagnetic  field.  To  describe  the  latter,  we  solve  the  Maxwell’s  equations  for 
propagation  of  electromagnetic  wave  (in  the  visible/near  IR)  through  the  plasma.  The  particle-in- 
cell  model  solves  the  relativistic  equations  of  motion  of  the  charged  particles  (ions  and  electrons) 
from  which  one  can  derive  detailed  information  about  the  particle  positions,  velocity  and  energy. 
Coupling  the  two  parts  of  the  model  turned  out  to  be  a  major  problem  since  the  plasma  density  is 
very  high.  Particle-in-cell  codes  tend  to  produce  “noisy”  quantities,  such  as  particle  densities  and 
currents.  The  numerical  noise  quickly  amplifies  and  the  solution  is  overwhelmed  by  artificial 
noise.  Therefore,  we  developed  a  technique  specifically  designed  to  couple  a  particle-in-cell  code 
to  a  Maxwell  solver.  We  tested  the  technique  for  both  underdense  and  overdense  plasma  and  for 
laser  intensities  in  the  entire  range  of  interest  (lO'®  -10^'  W/cm^).  The  tests  we  performed  were 
completely  satisfactory  in  terms  of  accuracy  and  robustness  and  indicate  that  the  technique  we 
developed  is  suitable  to  describe  high-intensity  laser-plasma  interactions  and  the  production  of 
high-energy  ion  beams  for  nuclear  fusion  and  neutron  production  (Ref  [2]). 

2. Ion  beam-target  deposition  model. 

A  Monte-Carlo  code  to  calculate  the  angular  distribution  of  neutrons  generated  by  a  beam 
of  deuterons  was  developed.  The  nuclear  fusion  reaction  considered  is  D  +  D  ^  n+^He  and  the 
target  is  deuterated  polyethylene  (CD2),  which  is  widely  used  in  this  kind  of  experiments.  The 
Monte-Carlo  code  takes  into  account  the  actual  trajectory  of  the  ion.  Ions  are  launched  with 
prescribed  initial  energy  and  direction  relative  to  the  target.  The  ion  loses  energy,  which  is 
accounted  for  by  the  energy  loss  function.  The  stopping  power  is  calculated  from  the  stopping 
powers  of  deuterons  in  carbon  and  deuterons  in  deuterium.  The  (high-energy)  deuterons  are  also 
scattered  due  to  collisions  with  the  target  material;  therefore  their  trajectories  are  not  a  straight 
line.  The  ion  energy  and  trajectory  are  followed  in  time  and  space  until  the  ion  energy  becomes 
less  than  a  prescribed  cut-off  energy  (~1  KeV).  During  each  time  step  the  neutron  flux  to  specific 
directions  of  observations  is  also  calculated  from  the  differential  cross  section.  The  model  outputs 
the  angular  distribution  of  neutrons,  as  well  as  the  total  neutron  yield. 
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3.  Coupling  of  the  laser-target  and  ion  beam-target  deposition  models. 

The  two  models  can,  in  general,  be  used  independently.  Here  they  are  connected  since  the 
output  of  the  laser-target  deposition  model  (deuteron  velocity  and  energy)  serves  as  input  of  the 
ion  beam-target  deposition  model.  The  laser-target  deposition  model  outputs  the  energy  and 
direction  of  every  deuteron  tracked  in  the  simulations,  which  is  precisely  the  input  for  the  Monte- 
Carlo  code.  The  output  of  the  Monte-Carlo  code  is  the  angular  distribution  of  neutrons  and  the 
total  neutron  yield  (in  absolute  units). 


Laser-target 

output 

tv 

deposition  model 

> 

Ei,  Gi 

(energy  and  angle 
of  each  ion) 


input 


J 


Ion  beam-target 
deposition  model 


Fig.  1  Coupling  of  the  laser-target  deposition  model  (the  electromagnetic  particle-in-cell  model)  to 

the  Monte  Carlo  ion  beam  target  deposition  model. 

B  Sample  results  and  accomplishments: 

1.  Double-layer  target 

The  formation  of  a  high-energy  deuteron  beam  from  the  interaction  of  ultra-short  laser 
pulses  with  a  planar  double-layer  solid  target  was  investigated.  The  target  is  a  two-layer  thin  foil 
with  thickness  of  ~1  pm.  The  front  layer,  which  accounts  for  most  of  the  target  mass,  is  made  of 
high  Z  material.  We  chose  gold,  which  is  commonly  used.  The  rear  layer  is  very  thin  (0.02  -  0.1 
pm)  and  it  is  made  of  deuterium-containing  material  (for  simplicity  we  assumed  that  it  is  made  of 
pure  deuterium).  Using  the  particle-in-cell  model,  we  studied  a  range  of  parameters  for  which  a 
highly  directional  deuteron  beam  is  created.  We  varied  the  following  laser-target  parameters:  (a) 
peak  laser  intensity:  from  10*^  to  10^'  W/cm^,  (b)  laser  pulse  duration:  from  40  to  160  fs,  and  (c) 
target  thickness:  from  0.1  and  3  pm.  The  goal  is  to  establish  trends  and  scaling  laws,  as  well  as  a 
suitable  range  of  (laser  and  target)  parameters  that  can  be  met  in  experiments.  Figure  2  summarizes 
two  important  quantities  related  to  laser  energy  deposition:  mean  energy  of  the  ion  beam  and 
conversion  efficiency  of  laser  energy  to  ion  energy.  The  regions  for  which  MeV  ions  are  generated 
are  shown  in  green.  Figure  3  plots  the  ion  beam  energy  and  angular  distribution  vs.  target 
thickness.  One  can  see,  for  example,  that  thinner  targets  produce  more  energetic  ions  (Fig.  3a,  top 
left).  Intermediate  thickness  favors  ion  beam  directionality  (Fig.  3d).  Large  target  thickness  (>  1 
pm)  leads  to  both  poor  directionality  (Fig.  3h)  and  low  energy  (Fig.  3g)  of  the  ion  beam. 
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Fig.  2.  Laser  energy  (top  row),  mean  energy  of  (seeond  row);  and  energy  conversion 
efficiency  of  Z)^  (bottom  row)  as  a  function  of  target  thickness  (left  column),  peak  laser  intensity 
(center  column)  and  laser  pulse  duration  (right  column). 
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Fig.  3  Deuteron  energy  distribution  function  (left  column)  and  angular  distribution  function  (right 
column)  for  different  target  thicknesses.  All  parameters  are  as  in  Fig.  2. 

The  next  step  is  to  calculate  the  neutron  yield  from  the  ion  beam-target  deposition  model.  The 
secondary  target  is  a  thick  (~  1mm)  slab  of  deuterated  material  (CD2).  The  Monte-Carlo  model 
follows  the  deuteron  trajectories  in  the  target  and  calculates  simultaneously  the  neutron  flux  to 
specific  directions  of  observations.  Figure  4  displays  a  typical  angular  distribution  function  of 
neutrons  from  a  CD2  target.  In  Fig.  4  we  plot  the  number  of  neutrons  per  unit  angle:  per  steradian 
(left)  and  per  radian  (right).  A  substantial  number  of  neutrons  are  emitted  at  small  angle  (0-20*^). 
The  total  number  of  neutrons  is  2x10^  at  laser  energy  of  1.1  J. 
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Fig.  4  Angular  distribution  function  of  neutrons  per  steradian  (left)  and  per  radian  (right). 
Figure  5  (below)  summarizes  our  parametrie  study. 
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Fig.  5  Average  deuteron  energy  (first  row);  neutron  yield  per  ion  (seeond  row),  total 
number  of  deuterons  (third  row)  and  total  neutron  yield  (last  row)  for  different  Au  layer 
thiekness  (first  column),  D  layer  thickness  (second  column),  peak  laser  intensity  (third 
column)  and  laser  pulse  duration  (right  eolumn).  Symbols  -  simulations  by  Toupin  et  al 
(Toupin  C,  Lefebvre  E,  and  Bonnaud  G  2001  Phys.  Plasmas  8  1011). 
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Figure  5  plots  the  average  incident  deuteron  energy  <  E‘"‘'  > ,  the  total  neutron  yield  per  ion,  the 

number  of  deuteron  ions  and  the  total  neutron  yield  are  plotted  vs.  thickness  of  the  gold  and 
deuterium  layers,  peak  laser  intensity  and  laser  pulse  duration.  For  the  intermediate  peak  laser 
intensity  of  1^  =  10^"'  W ! the  neutron  yield  is  10^  -10^  neutrons,  depending  on  the  laser  pulse 

duration.  The  results  shown  so  far  have  been  analyzed  in  detail  and  are  published  in  Plasma  Phys. 
Control.  Fusion  (Refs.  [3]). 

2.  Uniform  target 

The  target  design  of  the  double-layer  target  offers  a  lot  of  flexibility;  in  particular,  it  allows 
us  to  control  the  deuteron  energy  and  angular  distribution.  It  has,  however,  a  few  drawbacks. 
Complications  may  arise  due  to  alignment  of  the  two  targets,  etc.  Another  problem  is  the  disparity 
of  thicknesses  of  the  primary  and  secondary  targets  (primary  target  for  ion  acceleration,  secondary 
target  for  neutron  generation).  The  secondary  target  must  be  thick,  of  the  order  of  ~1  mm,  in  order 
to  be  able  to  stop  all  energetic  (~MeV)  deuterons,  while  the  primary  target  is  very  thin  (~1  pm), 
making  it  fragile  and  susceptible  to  damage  and  mis-alignment.  A  remedy  to  these  problems  is  to 
use  a  thicker  uniform  target,  which  can  be  made  of  the  same  material  as  the  secondary  target.  This 
approach  allows  consolidation  of  the  primary  and  secondary  targets  to  form  a  single  target.  The 
key  issue,  target  robustness,  is  now  achieved  by  the  target  thickness  (~1  mm).  This  is,  perhaps,  the 
main  reason  such  targets,  made  either  from  CD2  or  CD,  are  employed  in  virtually  all  experiments. 
The  objectives  are  to  study  the  ion  generation,  and  consequently,  the  neutron  generation  from  a 
single-layer  CD2  target  and  to  compare  it  with  the  neutron  generation  from  its  counterpart,  a 
double-layer  target.  This  will  help  us  design  better  targets  for  more  efficient  neutron  generation. 

Extensive  simulations  with  the  new  target:  a  planar  uniform  CD2  target  was  carried  out. 
Two  aspects  of  the  neutron  production  from  high-intensity  ultrashort  pulse  lasers  prompted  the 
present  study:  (i)  the  controversies  related  to  target  design  (material,  dimensions,  etc.);  and  (ii) 
studying  realistic  targets  employed  in  actual  experiments.  We  address  the  first  issue  by  comparing 
two  targets:  a  “thick”  (»1  pm)  uniform  piece  of  deuterated  plastic  (CD2)  and  a  double-layer 
“thin”  (<1  pm)  foil.  The  second  issue  is  addressed  by  comparing  simulation  results  with 
experimental  data.  First,  we  studied  the  deuteron  acceleration  mechanism  from  a  planar  uniform 
CD2  target.  It  originates  from  the  rear  target  surface  via  Target  Normal  Sheath  Acceleration  with  a 
modest  contribution  from  the  bulk  of  the  target. 
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Figure  6  displays  typical  energy 
and  angular  distribution 
functions  from  a  uniform  and 
double-layered  target  for  equal 
laser  parameters.  The  double¬ 
layered  target  produces  more 
energetic  ions  (a)  and  strongly 
forward-peaked  angular 
distribution  (c).  The  ion  energy 
and  angular  distributions  appear 
to  be  very  different  compared  to 
the  previous,  double  layer 
target. 
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Fig.  6  Energy  (a)  and  (b)  and  angular  (c)  and  (d)  distribution 
function  of  deuterons  from  a  double  layer  (Au+D)  target 
(left)  and  single  layer  (CD2 )  target  (right). 


However,  we  found  that  in  both  cases 
about  the  same  amount  of  energetic 
(above  1  MeV)  deuterons  are 
produced,  of  the  order  of  10'®-10'*. 
The  deuterons  from  the  double-layer 
target  (blue  triangles)  are  somewhat 
more  energetic  compared  to  the 
uniform  target  (red  open  squares). 
Figure  7  plots  the  maximum  deuteron 
energy  from  a  uniform  and  double¬ 
layer  targets. 
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Fig.  7  Maximum  deuteron  energy  vs.  laser  fluence  for 
double-layer  and  uniform  targets. 

The  main  difference  between  the  two  is  the  deuteron  angular  distribution:  thin  (sub-micron) 


double-layer  targets  produce  a  core  of  high-energy  collimated  (10-20°)  beam  of  deuterons,  while 
“thick”  (several  micron)  uniform  targets  produce  a  high-divergence  beam  of  deuterons  with 
somewhat  lower  energies. 
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Next,  we  investigated  the  neutron 


yield  in  absolute  units.  The  double-layer 
target  is  more  efficient,  but  with  the  laser 
energy  increasing  the  differences  between 
the  two  targets  diminish.  Above  about  10  J 
both  give  comparable  neutron  yields. 

Comparing  model  predictions  of  the 
neutron  yield  with  experimental  data 
benchmarked  the  simulation.  The  model 
predictions  are  in  good  agreement  with 
measured  neutron  yields  over  a  wide  range 
of  laser  energies,  spanning  two  orders  of 
magnitude. 


E 


laser 


(J) 


Fig.  8  Neutron  yield  from  a  uniform  CD2  target 
and  double  layer  target  vs.  laser  energy.  Solid  line 
with  symbols  -  model  predictions  for  double-layer 
and  uniform  targets;  symbols  -  experiments. 


Influence  of  the  preplasma 

3.  Target  design 

We  continue  our  investigations  focusing  on  the  target  design.  We  studied  double-layer  vs. 
uniform  targets  and  simultaneously  “thick”  ( » 1  pm)  vs.  “thin”  ( <  1  pm)  targets.  We  established 
that  double-layer  targets  tend  to  perform  better,  especially  in  terms  of  angular  distribution  of 
deuterons  and  neutrons.  All  previous  studies  have  been  limited  to  planar  targets  (Fig.  1  top),  which 
are  the  simplest  and  most  straightforward.  Shaping  targets  is  known  to  improve  the  ion  beam 
quality,  making  it,  for  example,  more  laminar  or  mono-energetic.  Therefore  our  next  goal  is  to 
perform  simulations  with  more  complex  (in  terms  of  shape)  targets. 

We  will  first  study  the 
impact  of  the  pre-plasma  on  the 
neutron  production.  The  pre-plasma 
occurs  during  the  interaction  of  a 
planar  target  with  the  foot  of  the 
laser  pulse.  This  laser  pre-pulse  is 
always  present  and  may  change  the 
target  geometry  by  ionizing  the 
irradiated  front  surface  of  the  target. 

A  layer  of  plasma  is  formed,  called 
pre-plasma,  the  extent  of  which 
depends  upon  the  parameters  of  the 
laser  pre-pulse.  Figure  1  top  shows 
an  “ideal”  planar  target,  unaffected 
by  the  laser  pre-pulse.  Figure  1 
bottom  shows  the  more  realistic 
situation  of  a  planar  target  with  a 
pre-plasma.  The  pre-plasma  length 
can  be  (at  least  conceptually) 
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controlled  to  some  degree  by  Fig.  1  Planar  target  (top)  and  a  target  with  a  pre-plasma 
varying  the  intensity  of  the  laser  pre-  (bottom).  The  thickness  of  the  CD2  layer  is  50  nm. 
pulse. 


The  neutron  yield  in  absolute 
units  vs.  pre-plasma  length  is  shown  in 
Fig.  2  for  three  thicknesses  of  the  Au 
layer.  The  points  to  the  left  at  0.02  pm 
mimic  “ideal”  planar  target.  As  it  can  be  ^ 
seen  there  is  an  optimum  pre-plasma  c 
thickness  of  the  order  of  0. 1-0.3  pm,  for 
which  the  neutron  yield  peaks.  So  the 
presence  of  some  pre-plasma  is  actually 
helpful  for  the  neutron  production. 

The  thicker  the  target,  the  more 
pronounced  is  the  pre-plasma  impact 
(green  line  &  triangles,  corresponding  to 
Au  layer  of  1pm).  The  neutron  yield  can 
increase  by  a  factor  of  ~  100  compared  to 
“ideal”  targets  with  no  pre-plasma.  In  any 
case  the  pre-plasma  (and  the  target  shape 
in  general)  appears  to  be  important  for 
the  neutron  production. 


Fig.  2  Neutron  yields  from  a  double-layer  AU-CD2  target 
vs.  pre-plasma  length.  The  three  curves  are  for  a  thickness 
of  the  Au  layer  of  0.25  pm  (red  square),  0.50  pm  (blue 
cycles),  and  1  pm  (green  triangles). 
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Abstract 

The  formation  of  a  high-energy  ion  beam  suitable  for  driving  nuclear  fusion 
reactions  and  producing  MeV  neutrons  was  investigated.  The  interaction  of 
intense  ultra  short  laser  pulses  with  a  double-layer  thin  foil  for  the  production 
of  MeV  deuterons  was  studied  theoretically  and  numerically  simulated  using 
a  two-dimensional  electromagnetic  particle-in-cell  model.  The  directionality 
and  energy  of  the  deuteron  beam,  specifically  the  conversion  efficiency  of 
laser  energy  into  deuteron  kinetic  energy,  and  deuteron  energy  and  angular 
distribution  functions  are  studied  as  a  function  of  peak  laser  intensity,  laser 
pulse  duration  and  target  thickness.  A  range  of  parameters  was  determined  for 
which  a  highly  directional  deuteron  beam  is  generated. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

High-intensity  ultrashort  pulse  lasers  are  compact  and  versatile  systems  for  inducing  nuclear 
fusion  reactions.  It  has  been  known  for  some  time  that  when  high-intensity  laser  radiation 
interacts  with  solid  targets,  it  creates  energetic  particles,  capable  of  driving  fusion  reactions. 
The  remarkable  progress  in  the  development  of  laser-based  particle  accelerators  is  pushing 
the  frontiers  in  basic  research  and  holds  promise  for  harnessing  nuclear  fusion  for  practical 
applications.  As  laser  systems  producing  ultrashort  high-intensity  pulses  are  expected  to 
become  more  compact  and  affordable  several  advanced  applications  can  be  envisioned,  such  as 
positron  emission  tomography  [1],  cancer  therapy  [2,3]  and  laser-induced  nuclear  reactions  [4]. 
Laser-driven  neutron  sources  are  an  alternative  to  the  accelerator-  and  reactor-driven  sources 
offering  high  brightness,  compactness,  short  duration,  shielding  and  relatively  low  cost.  This  is 
particularly  advantageous  for  applications  such  as  fast  neutron  radiography  [5],  transmutation 
of  nuclear  waste  [6]  and  fusion  research  [7]. 
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Figure  1.  Neutron  production  concept  (a)  and  target  configuration  (b).  The  laser  pulse  propagates 
from  left  to  right  along  the  axis  ‘x’,  the  laser  electric  and  magnetic  fields  are  along  the  ‘y’  and  ‘z’ 
axes,  respectively. 


For  our  purposes  neutron  production  can  be  visualized  schematically  as  shown  in 
figure  1(a).  During  the  first  stage,  a  high-energy  ion  beam  is  generated.  In  the  next  stage, 
fast  neutrons  are  born  in  fusion  reactions  driven  by  the  fast  ions.  During  the  third  stage,  the 
neutrons  impinge  on  the  target  of  interest  and  the  interaction  of  the  fast  neutrons  with  the  target 
leads  to  the  emission  of  gamma  rays.  The  gamma-ray  ‘signature’  is  unique  for  each  element 
and,  in  principle,  can  be  used  to  identify  the  constituents  of  the  material  producing  them.  The 
real  challenge  is  to  produce  a  neutron  beam  in  sufficient  numbers  (10^-10^^)  that  is  reasonably 
well  collimated  (10°-20°  divergence).  The  difficulty  is  rooted  in  the  precursor,  the  ion  beam, 
which  must  be  both  energetic  and  highly  directional.  The  ion  energy  must  be  comparable  or 
higher  than  the  knock-off  neutron  energy  (typically  "^2.5  or  "^14  MeV).  Only  then  the  neutrons 
bom  in  nuclear  fusion  reactions  are  scattered  preferentially  in  the  forward  direction.  Thus  we 
focus  our  attention  on  laser-target  interactions  capable  of  producing  energetic  deuteron  beams 
with  small  divergence. 

Ion  acceleration  can  be  achieved  by  various  types  of  laser-target  interactions:  underdense 
plasmas,  clusters  and  solids.  Generation  of  multi-MeV  ion  beams  with  narrow  angular 
spread  has  been  demonstrated  in  the  interaction  of  high-intensity  ultrashort  pulse  lasers  with 
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underdense  plasma  [8, 9].  Neutron  production  from  nuclear  fusion  reactions  in  laser-cluster 
interactions  has  also  been  demonstrated  [10].  The  most  attractive  feature  of  this  approach  is 
the  use  of  off-the-shelf  laser  systems  as  nuclear  fusion  reactions  can  be  triggered  at  peak  laser 
intensities  as  low  as  "^10^^  W  m~^  ("^10^^  W  cm“^).  Unfortunately,  the  neutrons  created  by 
laser-cluster  interactions  are  scattered  isotropically  [11]  and  their  number  is  only  of  the  order 
of  10"^- 10^  per  shot  [12, 13].  Therefore  clusters  are  not  suitable  for  generation  of  an  energetic 
focused  beam  of  neutrons.  Laser-driven  ion  acceleration  and  generation  of  fast  neutrons  can 
also  be  achieved  as  a  result  of  ultraintense  laser  irradiation  of  solid  targets  [14,  15].  It  is 
now  well  established  that  ions  at  the  rear  surface  of  a  thin  (~  1  /xm)  foil  can  be  accelerated  to 
relativistic  energies  by  the  electrostatic  field  of  a  space-charge  potential  created  by  escaping 
non-thermal  electrons  and  the  ions  that  are  left  behind.  This  mechanism  is  known  as  target 
normal  sheath  acceleration  (TNSA)  [15-17].  Recent  work  by  Sentoku  et  al  elaborated  on 
this  widely  accepted  mechanism  of  ion  acceleration  from  the  rear  surface  of  the  target.  They 
found  that  the  effect  of  electron  recirculation  enhances  the  peak  proton  energy  if  the  target 
thickness  is  shorter  than  half  of  the  laser  pulse  length  [18].  When  the  target  thickness  becomes 
comparable  to  the  plasma  skin  depth,  the  laser  penetrates  to  the  rear  of  the  target  and  TNSA 
transitions  to  the  so-called  ‘laser  breakout  afterburner’  regime  [19].  This  regime  distinguishes 
three  stages  of  ion  acceleration:  a  period  of  TNSA,  followed  by  a  period  of  enhanced  TNSA 
during  which  the  cold  electrons  are  converted  into  hot  electrons  and  rapid  ion  acceleration  in  the 
enhanced  longitudinal  electrostatic  field.  At  intensities  exceeding  10^^  W  m~^  laser-induced 
shocks  could  also  accelerate  ions  to  high  energies  [20] .  At  even  higher  intensities,  for  which  the 
dimensionless  laser  field  amplitude  becomes  comparable  to  the  ratio  of  ion-to-electron  mass, 
direct  laser  acceleration  of  protons  to  relativistic  energies  becomes  possible  [21].  Regardless 
of  the  acceleration  mechanism,  if  the  foil  parameters  are  properly  chosen,  a  well-collimated 
beam  of  light  ions  can  be  accelerated  to  energies  in  the  MeV  range.  Recent  progress  in 
ultrashort  pulse  laser  technology  now  allows  generation  of  laser  pulses  with  intensities  which 
are  adequate  for  producing  an  ion  beam  with  the  required  parameters.  In  spite  of  the  practical 
challenges  of  implementing  such  laser  systems,  both  the  ion  beam  collimation  and  high  energy 
are  crucial  for  production  of  high- flux  neutron  beam. 

Neutron  generation  in  high-intensity  laser-foil  interaction  has  been  demonstrated 
experimentally  [22-31]  and  theoretically  [7, 32-35].  Perkins  etal  [7]  surmised  various  regimes 
of  laser-target  interactions  depending  on  target  configuration  and  peak  laser  intensity.  Of 
particular  interest  are  the  results  by  Toupin  et  al  [33],  which  studied  the  angular  distribution 
of  neutron  emission  as  a  function  of  laser  intensity  and  maximum  electron  density.  Their 
results  indicate  that  if  the  medium  generating  light  ions  is  highly  overdense  plasma  (typical 
for  laser-thin  foil  interaction),  the  neutrons  are  emitted  preferentially  in  the  forward  direction. 
Macchi  et  al  [35]  studied  the  acceleration  of  ion  bunches  by  using  an  analytical  model  as  well 
as  one-  and  two-dimensional  particle-in-cell  (PIC)  models.  They  were  concerned  mostly  with 
the  formation  of  high-density  low-energy  (<  1  MeV)  ion  beams,  which  are  not  optimal  for  the 
formation  of  directional  neutron  beams.  To  our  knowledge,  there  is  no  systematic  theoretical 
study  regarding  the  formation  of  collimated  deuteron  beams,  except  for  [33-35] .  Our  objective 
is  to  study  theoretically  the  formation  of  a  collimated  beam  of  high-energy  (MeV)  deuterium 
ions  as  a  result  of  the  interaction  of  high-intensity  laser  radiation  with  thin  double-layer  foil. 
This  is  considered  as  the  first  stage  of  a  general  scheme  leading  to  the  generation  of  directional 
neutron  beams.  Our  primary  interest  is  in  studying  the  directionality  and  energy  of  the  deuteron 
beam,  specifically  (i)  the  conversion  efficiency  of  laser  energy  into  ion  kinetic  energy,  (ii)  ion 
energy  distribution  function  and  (iii)  ion  angular  distribution  function.  These  data  will  later 
serve  as  input  parameters  of  a  model  simulating  the  production  and  scattering  of  the  neutron 
beam.  In  section  2  we  outline  the  target  configuration  and  the  approach  we  use  to  model  the 
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laser-target  interaction.  In  section  3  we  report  on  results  for  a  specific  target  and  summarize 
our  findings  in  section  4. 

2.  Laser-target  interaction  model 

To  simulate  the  interaction  between  intense  laser  radiation  and  solid  target  we  implemented  a 
2D  relativistic  electromagnetic  PIC  method.  The  target  is  approximated  as  a  slab  of  overdense 
plasma  whose  constituents  (electrons  and  ions)  are  modeled  with  the  PIC  method.  The 
propagation  of  laser  radiation  into  the  plasma  is  modeled  by  solving  the  Maxwell  equations. 
Other  numerical  models  that  adopted  this  approach  have  been  developed  for  2D  [35-56]  and 
3D  [21, 57-62]  geometry.  Our  version  of  a  2D  model  is  given  in  appendix,  where  the  coupling 
of  the  electromagnetic  solver  to  the  PIC  code  is  described.  Here  we  will  provide  only  a  brief 
account  for  the  reasons  that  prompted  us  to  develop  this  algorithm.  The  most  obvious  approach 
to  couple  the  PIC  method  to  a  Maxwell  solver  is  to  calculate  the  current  density  from  PIC, 
substitute  it  into  the  Maxwell  equations  and  solve  them  for  the  electromagnetic  field.  This  two- 
step  procedure  is  repeated  for  each  time  step.  This  intuitive  approach,  however,  suffers  from 
some  deficiencies.  PIC  codes  are  notorious  for  producing  ‘noisy’  quantities,  such  as  particle 
densities  and  currents.  The  numerical  noise  quickly  amplifies  and  spurious  numerical  solutions 
are  generated.  Therefore  strategies  that  exploit  the  specific  nature  of  the  PIC  method  must  be 
sought  to  remedy  the  problem  of  numerical  instability.  We  implemented  a  relationship  between 
the  laser  electric  field  E  and  current  density  j  generated  by  the  movement  of  charged  particles 
for  the  special  case  when  j  is  calculated  from  a  2D  PIC  model  in  Cartesian  coordinates: 

jn+\/2  ^  -«+l/2^n+l/2 

Equation  (1)  is  concisely  written  in  matrix  form,  in  which  a  is  a  2  x  2  matrix  with  elements 
calculated  from  quantities  such  as  particle  mass,  charge,  density,  relativistic  factor  and 
cyclotron  frequency  (appendix).  The  superscript  denotes  the  discretization  time  level  with 
'n'  being  the  previous  and  n  1 — the  new  time  level.  If  the  current  density  (1)  with  the 
electric  field  factored  out  is  substituted  into  the  Maxwell  equations,  the  numerical  scheme 
becomes  very  robust.  Numerical  noise  and  grid  heating  are  strongly  suppressed  and  spurious 
numerical  solutions  are  eliminated.  The  technique  is  applicable  to  modeling  the  interaction  of 
laser  radiation  with  both  underdense  and  overdense  plasmas,  but  it  is  particularly  useful  for 
highly  overdense  plasma  (~100  times  the  critical  electron  density),  where  other  techniques 
have  difficulties  or  fail.  Decomposition  of  the  current  density  in  the  form  (1)  has  also  been 
used  in  other  electromagnetic  PIC  models,  but  the  specific  numerical  implementation  may 
have  been  different  from  ours  [63, 64]. 

The  target  is  a  two-layer  thin  foil  with  thickness  between  0.1  and  3  /xm.  The  front  layer, 
which  accounts  for  most  of  the  target  mass,  is  made  of  high-Z  material.  The  rear  layer 
is  very  thin  (0.01-0.1  /xm)  and  it  is  made  of  deuterium-containing  material.  The  choice  of 
high-Z-low-Z  material  is  motivated  by  the  observation  that  if  the  ratio  of  mass  to  charge 
for  the  front  and  rear  material  is  sufficiently  large,  the  light  ions  from  the  rear  surface  are 
accelerated  much  more  efficiently  than  the  heavy  ones  [57].  For  definitiveness  we  chose  gold 
as  a  front  layer  material,  which  is  a  common  substance  in  this  kind  of  high-intensity  laser- 
target  interactions.  The  plasma  slab  corresponding  to  this  layer  has  a  thickness  Lau,  density 
fXAu  and  charge  Zau  —  L  respectively.  For  simplicity  we  assumed  that  the  rear  layer  is  made 
of  pure  deuterium  with  thickness  and  density  Ld  and  fxo,  respectively.  Both  the  front  and  the 
rear  layers  have  the  form  of  a  disc  with  a  prescribed  diameter  Dau  and  Dd,  respectively. 
The  geometry  of  the  target  is  shown  in  figure  \{h).  More  complex  targets  in  terms  of 
composition  and  geometry  will  be  modeled  in  the  future.  The  laser  radiation  is  an  ultrashort 
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pulse  of  linearly  polarized  light  normally  incident  on  the  target.  The  laser  pulse  has  a 
spatio-temporal  prohle  I{y,t)  =  /q  sin^(7r^/To)  exp(— (y/7?o)^)  with  /q  being  the  peak  laser 
intensity,  Tq  the  laser  pulse  duration  and  Rq  the  laser  spot  radius.  The  laser  propagation 
direction  is  -Fx.  The  simulation  box  is  a  rectangle  with  length  (along  the  laser  propagation 
direction)  and  width  Ly  (in  the  transverse  direction).  There  is  a  vacuum  region  in  front  of 
the  target  with  a  length  xau  extending  up  to  several  laser  wavelengths.  A  vacuum  region  of 
comparable  length  behind  the  target  is  also  present.  The  computational  domain  is  partitioned 
into  a  suitable  number  of  grid  cells  x  Ny  with  cell  size  Ax  x  Ay,  where  Ax  =  L^/Nx 
and  Ay  =  Ly/Ny.  The  cell  size  is  of  critical  importance.  On  one  hand,  it  cannot  be  made  too 
small  since  for  a  given  number  of  particles  there  would  be  too  few  particles  per  cell,  which 
results  in  poor  statistics.  On  the  other  hand,  too  large  a  cell  size  may  lead  to  an  incorrect 
solution  of  the  Maxwell  equations  since  the  skin  depth  must  be  resolved.  Therefore,  Ax  is 
chosen  to  be  0. 1-0.5  times  the  skin  depth  and  in  most  calculations  it  is  10  nm.  The  number  of 
particles  Ap,  used  in  the  PIC  model,  is  typically  a  few  million.  The  time  step  At  was  chosen 
based  on  the  following  consideration.  We  found  empirically  that  the  best  results  are  obtained 
if  a  particle  resides  in  a  given  cell  for  several  time  steps.  This  condition  may  be  expressed  as 
At  ^  (1/2) Ax /max {\Yk\},  where  the  denominator  is  twice  the  expected  maximum  velocity 
of  any  particle.  In  the  highly  relativistic  limit  this  condition  reduces  to  At  ^  Ax /2c,  where  c 
is  the  speed  of  light. 

3.  Results  and  discussions 

Electromagnetic  PIC  simulations  are  performed  for  a  range  of  laser  and  target  parameters 
dehned  in  the  previous  section.  As  mentioned  in  the  introduction,  we  are  primarily  interested 
in  the  conversion  efficiency  of  laser  energy  into  ion  kinetic  energy  and  the  ‘quality’  of  the 
ion  beam  (energy  and  angular  distribution  functions),  but  other  related  parameters  such  as 
mean  electron  energy  and  electron  energy  distribution  function  will  also  be  discussed.  We 
will  commence  our  analysis  with  the  temporal  characteristics  of  a  typical  case  (section  3.2). 
In  section  3.3  we  will  vary  laser  and  target  parameters  in  order  to  derive  scaling  laws  and 
parameter  trends. 

3.1.  PIC  simulation  parameters 

The  peak  laser  intensity  and  wavelength  for  the  typical  case  are  /q  =  lO^^^Wm"^  and 
ko  =  1  /xm,  respectively.  The  transverse  laser  prohle  is  a  Gaussian  with  l/e  spot  size 
/?o  =  3  /xm.  The  laser  pulse  duration  is  Tq  =  80  fs  with  FWHM  Tfwhm  =  ro/2  of  40fs  (12 
cycles)  and  the  input  laser  energy  £^iaser  =  fo^  f)dt'  =  tt/^q/oTfwhm  is  1.13  J.  The 
computation  time  is  extended  to  160  fs,  twice  the  laser  pulse  duration,  since  the  deuterium  ions 
continue  their  acceleration  after  the  end  of  the  laser  pulse.  The  target  is  a  highly  overdense 
plasma  with  density  =  lOOfZc  (n^  =  1.12  x  10^^  m“^ — critical  electron  density).  The 
gold  and  deuterium  layers  have  thicknesses  of  0.25  and  0.05  /xm,  respectively.  The  width 
of  both  layers  is  twice  the  laser  spot  diameter,  e.g.  Dau  =  =  12  /xm.  As  the  laser 

intensity  drops  signihcantly  at  the  periphery  of  the  target  (~l/50  of  that  on  the  axis),  the 
contribution  of  particles  located  at  y  >  2/?o  to  the  energy  absorption  and  ion  beam  acceleration 
is  small.  Doubling  the  target  width  was  found  to  have  a  negligible  impact  on  both  the  energy 
absorption  and  angular  distribution  of  deuterons.  Detailed  information  about  the  laser  and 
target  parameters,  as  well  as  some  additional  simulation  parameters,  is  given  in  table  1.  The 
laser  pulse  prohle  /(y  =  0,  t)  and  the  unperturbed  laser  electric  held  =  0,  t)  cos((X>oO 
on  axis  with  Eq  =  ^(21  /s^c)  being  the  envelope  of  the  laser  electric  held  and  coq  the  central 
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Table  1.  Laser,  target  and  computational  parameters. 


Parameter 

Value 

Laser  peak  intensity  Iq 

Laser  relativistic  parameter  ao 

8.5 

Laser  wavelength  Xq 

1  /xm 

Laser  period  Tq  =  Xq  jc 

3.33  fs 

Laser  focal  spot  radius  Rq 

3  /xm 

Laser  FWHM  tfwhm 

40  fs 

Laser  pulse  duration  to 

80  fs 

Front  layer  Lau  x  Dau 

0.25  X  12  /xm^ 

Front  layer  density  ^au 

5.8  X  10^8  m-3 

Front  layer  charge  Zau  —  1 

2 

Rear  layer  Ld  x  Dd 

0.05  X  12 /xm^ 

Rear  layer  density  no 

4.8  X  10^^  m-3 

Rear  layer  charge  Zd  —  1 

1 

Computation  time  T 

160  fs 

Simulation  box  Lx  Ly 

10  X  20  ixw3 

Grid  cells  Nx  x  Ny 

1000  X  2000 

Cell  size  Ax  x  Ay 

10  X  10  nm^ 

Skin  depth 

16  nm 

Length  of  vacuum  region  xau 

4  /xm 

Number  of  particles 

2  X  10'’ 

Time  step  Ar 

5  X  10-3  To 

cAt/Ax 

0.5 

laser  frequency  are  given  in  figure  2.  At  time  ^  =  0  the  front  of  the  laser  pulse  is  located  at 
spatial  position  v  =  xau,  the  vacuum-front  layer  interface. 

3.2.  Time -dependent  results 

In  our  first  example  we  will  focus  on  the  temporal  characteristics  of  the  electrons  and  the 
ions.  One  of  the  main  issues  of  the  laser-target  interaction  process  is  the  absorbed  laser 
energy  by  the  target.  When  a  laser  pulse  propagates  through  a  solid  medium  there  is  an 
energy  loss  associated  with  coupling  of  the  laser  electromagnetic  wave  to  the  plasma.  The 
absorption  of  laser  energy  is  plotted  in  figure  2(c).  The  solid  line  is  the  amount  of  input  laser 
energy  £^iaser(0  =  /o  ^  ’  which  enters  the  target  from  the  onset  of  the  laser  pulse 

t'  =  0  up  to  time  t'  =  t.  The  other  three  curves  are  the  energy  absorbed  by  electrons,  gold 
and  deuterium  ions.  They  comprise  the  energy  contribution  of  all  particles  belonging  to  the 
corresponding  specie,  regardless  of  their  position  and  direction  of  motion.  The  energy  absorbed 
by  electrons  peaks  near  the  end  of  the  laser  pulse  and  then  slightly  decreases.  The  decrease 
is  associated  with  a  transfer  of  kinetic  energy  from  the  electrons  to  the  ions.  An  additional 
(unwanted)  source  of  energy  loss  is  the  escape  of  particles  from  the  computational  box.  At  the 
end  of  the  simulations  the  electrons  have  absorbed  about  ~0. 1 1  J,  gold  ions — "^0.15  J  and  the 
deuterons — ~0.07  J.  The  total  absorbed  energy  of  "^0.33  J  represents  ~30%  of  the  input  laser 
energy.  The  conversion  efficiency  of  laser  energy  into  deuteron  kinetic  energy  is  ~6%.  The 
magnitude  of  the  energy  absorbed  by  the  front  and  rear  layers,  0. 15  J  and  ~0.07  J  respectively, 
deserves  further  elaboration.  The  front  layer  is  five  times  thicker  compared  with  the  rear  layer 
and  contains  more  atoms.  In  addition,  the  laser  field  decays  on  a  distance  of  a  few  skin  depths 
barely  reaching  the  rear  layer,  therefore  depositing  electromagnetic  energy  preferentially  in 
the  front  layer.  But  the  rear  layer  has  the  advantage  of  having  light  ions,  which  are  more  likely 
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Figure  2.  Laser  pulse  intensity  (a)  and  laser  electric  field  (b)  versus  time.  Absorbed  energy  by 
electrons  (dashed  line),  Au^"^  (open  cycles)  and  (open  squares)  and  input  laser  energy  (solid  line) 
(c)  versus  time.  The  laser  and  target  parameters  are:  Iq  =  10^^  W  m~^,  A.o  =  1  /xm,  Rq  =  3  /xm, 
tq  =  80 fs,  Fiaser  =  1.13  J,  Lau  =  0.25 /xm,  Ld  =  0.05 /xm,  =  lOOric.  The  simulation 
parameters  are  listed  in  table  1 . 


to  be  accelerated  than  the  heavy  gold  ions  from  the  front  layer.  At  these  conditions  both  layers 
absorb  about  equal  amounts  of  energy. 

Details  of  the  acceleration  process  are  given  in  figure  3,  in  which  we  make  a  direct 
comparison  of  the  properties  of  each  species.  The  top  row  plots  the  average  energy  of  the  species 
under  consideration.  The  peak  average  kinetic  energies  are  comparable,  of  the  order  of  1  MeV. 
The  corresponding  energy  distribution  functions  are  plotted  in  the  next  row  of  figures.  For 
both  electrons  and  ions  the  energy  distribution  function  is  non-Max wellian,  exhibiting  a  long 
high-energy  tail.  The  third  row  of  figures  plots  the  average  velocities  (Vq;)  =  (V^q;)  ^^,0; 

of  specie  «  =  e,  Au^"^,  along  the  laser  propagation  direction  v  and  the  transverse  direction 
y.  The  sum  is  over  all  particles  k  belonging  to  specie  a  and  Na  is  the  total  number  of  particles 
of  specie  a.  Figure  3(g)  clearly  shows  strong  oscillations  of  the  electrons  in  the  transverse 
direction,  parallel  to  the  laser  electric  field,  with  frequency  equal  to  that  of  the  fundamental 
laser  frequency.  A  small  positive  push  in  the  laser  propagation  direction  is  also  observed.  The 
mean  velocity  {Vx,e)  oscillates  with  frequency  equal  to  twice  the  fundamental  laser  frequency 
consistent  with  the  y  x  B  part  of  the  Lorentz  force.  The  longitudinal  ion  average  velocities 
look  quite  different.  The  magnitude  of  the  average  velocity  of  Au^"^  is  only  10^  m  s“^  and  the 
average  distance  the  gold  ions  can  move  during  the  short  laser  pulse  duration  of  "^10^  fs  is  only 
0.01  /xm.  Therefore,  for  all  practical  purposes  the  gold  ions  can  be  considered  immobile  on 
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Figure  3.  Mean  energy  of  electrons  (a),  (b)  and  (c);  energy  distribution  function  of 

electrons  (d),  Au^"^  (e)  and  (/);  mean  velocity  of  electrons  (g),  Au^"^  (h)  and  (/);  and  angular 
distribution  function  of  electrons  (j),  Au^+  (k)  and  D+  (/).  All  parameters  are  as  in  figure  2. 


such  a  short  time  scale.  It  should  be  kept  in  mind,  however,  that  they  absorb  a  non-negligible 
amount  of  laser  energy.  The  most  relevant  information  regarding  the  average  velocities  of  ions 
is  illustrated  in  figure  3(/).  The  deuterons  are  accelerated  to  an  average  velocity  approaching 
10^ms“\  which  is  {Vx,d+)/c  =  0.025  relative  to  the  speed  of  light.  The  angular  scattering 
is  expressed  quantitatively  through  the  angular  scattering  distribution  function  O(^),  where 
the  angle  (p  is  accounted  for  from  the  direction  of  laser  propagation  (figures  3(/)-(0)-  The 
angular  distribution  function  of  specie  a  is  constructed  by  calculating  the  angle  cpk  for  each 
particle  k  belonging  to  specie  a  according  to  cpk  =  \yy,k\/^x,k  and  then  making  a  histogram  of 
the  number  of  counts  versus  angle  bins.  The  histograms  in  figure  3  are  normalized  to  unity, 
i.e.  /q  =  1.  As  seen  in  figure  3(/)  the  electrons  are  uniformly  scattered,  while  the 

angular  spread  of  virtually  all  deuterons  lies  within  a  cone  of  about  20°  (figure  3(/)).  The 
angular  scattering  of  may  be  explained  from  geometrical  considerations.  The  problem 
is  essentially  one-dimensional  since  the  aspect  ratios  for  both  gold  and  deuterium  layers  are 
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Figure  4.  Snapshots  of  particle  positions  of  electrons  (left  column),  Au^'*'  (center  column)  and  D'*' 
(right  column).  All  parameters  are  as  in  figure  2. 


much  smaller  than  unity  (Lau/^Au  =  2  x  10“^  and  Ld/Z^d  =  4  x  10“^,  respectively).  It  is 
reasonable  to  expect  that  the  ions  will  be  accelerated  in  a  direction  normal  to  the  surface.  On 
the  other  hand,  transverse  inhomogeneities  of  the  accelerating  electric  field  cause  some  small 
angular  spread  of  the  accelerated  deuterons. 

The  main  features  of  electron  and  ion  acceleration  are  visualized  in  figure  4,  where 
snapshots  of  the  vy  positions  of  the  particles  are  shown.  The  left  column  refers  to  electrons, 
the  middle — to  gold  ions  and  the  right — to  deuterium  ions.  The  first  row  plots  the  particle 
positions  at  the  beginning  of  the  simulations.  The  second  row  refers  to  time  t  =  40  fs. 
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which  corresponds  to  the  peak  of  the  laser  pulse  intensity.  At  this  time  no  significant  particle 
movement  is  observed,  though  there  are  a  few  scattered  electrons.  The  next  row  of  figures 
shows  the  particle  positions  at  the  end  of  the  laser  pulse  {t  =  80  fs).  We  observe  electrons 
scattered  both  in  forward  and  backward  directions.  The  backward  scattered  electrons  are 
ejected  from  the  front  gold  surface.  The  forward  accelerated  electrons  are  a  result  of  the 
\  X  B  component  of  the  Lorentz  force  pushing  the  energetic  electrons  in  the  direction  of  laser 
propagation.  It  is  most  effective  on  axis  (y  =  0),  where  the  magnetic  field  is  strongest.  The 
electrons  accelerated  forward  are  relativistic  and  can  reach  the  rear  side  of  the  target.  They 
become  separated  from  the  immobile  gold  ions  inducing  a  space-charge  electrostatic  field, 
which  accelerates  the  deuterons.  The  most  energetic  deuterons  have  gained  momentum  and 
advanced  forward  a  distance  of  1  /xm.  The  snapshots  of  the  particle  positions  at  later  time 
t  =  120fs  and  at  the  end  of  the  simulations  t  =  160fs  show  that  there  is  a  cloud  of  blow- 
off  energetic  electrons  separated  from  the  target  and  accelerated  forward,  as  well  as  a  cloud 
of  deuterons  following  these  electrons.  The  latter  exhibits  a  bell-like  profile,  also  observed 
in  [35].  This  profile  is  a  result  of  the  radial  dependence  of  the  accelerating  electrostatic  field, 
which  is  largest  on  axis.  The  electrons  (and  ions)  are  preferentially  accelerated  on  the  axis  and 
go  ahead  of  those  at  the  periphery,  forming  the  bell-like  shape  seen  in  figure  4. 


3.3.  Variation  of  the  peak  laser  intensity,  laser  pulse  duration  and  target  thickness 

In  this  section,  we  describe  the  general  features  of  deuteron  acceleration  as  a  function  of  peak 
laser  intensity,  laser  pulse  duration  and  thickness  of  the  front  and  rear  layers.  Parametric 
studies  have  been  performed  by  others  [15,  18,48,65-67].  Our  analysis  is  limited  to  the 
three  most  important  parameters  related  to  energy  absorption  by  the  target:  absorbed  energy 
and  its  distribution  among  various  degrees  of  freedom,  deuteron  mean  energy  and  energy 
conversion  efficiency.  The  deuteron  acceleration  is  then  analyzed  based  on  the  energy  and 
angular  distribution  functions.  When  one  of  the  input  parameters  is  varied  (for  example,  target 
thickness),  all  other  simulation  parameters  are  kept  the  same.  By  ‘target  thickness’  we  mean 
the  extent  of  the  front  layer.  The  thickness  of  the  rear  layer,  usually  tens  of  nanometers,  is 
typically  much  smaller  compared  with  that  of  the  front  layer.  This  is  almost  always  the  case 
because  in  actual  experiments  the  light  ions  (H"^)  are  just  contaminants  residing  on  the  back 
surface  of  the  target. 

We  commence  with  variation  of  the  target  thickness.  To  understand  the  incentives  of 
studying  the  variation  of  this  parameter,  let  us  compare  two  commonly  used  targets  in  laser- 
target  interactions:  clusters  and  solids.  Clusters  are  tiny  targets  at  solid  density,  only  a  few 
nanometers  in  diameter.  Extensive  studies  during  the  last  ten  years  revealed  that  the  ion 
kinetic  energy  of  exploding  clusters  increases  dramatically  with  cluster  diameter  (which,  for 
simplicity,  may  be  attributed  to  as  ‘target  thickness’)  [68].  On  the  other  hand,  when  the  target 
becomes  a  few  micrometer  thick,  the  energy  of  the  light  ions  decreases  with  target  thickness 
increasing  as  shown  by  MacKinnon  et  al[Al].  Apparently,  the  peak  ion  energy  is  expected  to 
reach  a  maximum  somewhere  in  between.  The  first  series  of  simulations,  in  which  the  target 
thickness  is  varied  between  0.1  and  3  /xm,  is  plotted  in  figure  5,  left  column.  We  observe 
two  distinct  cases  of  target  thickness.  Very  thin  front  layer  Lau  ^  0.5  /xm  provides  efficient 
energy  absorption  (figure  5 (a)),  high  deuteron  mean  ion  energy  ~lMeV  (figure  5(^))  and 
significant  energy  conversion  efficiency  of  a  few  percent  (figure  5(/)).  The  energy  conversion 
efficiency  y/D  is  defined  as  the  energy  absorbed  by  divided  by  the  input  laser  energy. 
Targets  with  thickness  Lau  ^  0.5  /xm  lead  to  a  precipitous  drop  of  both  the  mean  energy  and 
energy  conversion  efficiency  of  D"^.  The  reason  becomes  apparent  when  examining  the  specie 
energy  absorption  plotted  in  figure  5{a).  For  the  thickest  target  considered  (Lau  =  3  /xm)  the 
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Figure  5.  Laser  and  absorbed  energy  ia)-{d)\  mean  energy  (e)-(h)  and  energy  conversion  efficiency 
of  D+  (0-(0  as  a  function  of  target  thickness  (left  column),  D  layer  thickness  (second  column), 
peak  laser  intensity  (third  column)  and  laser  pulse  duration  (right  column). 


bulk  of  the  absorbed  energy  (~99%)  goes  to  the  ions  of  the  front  layer,  while  the  deuterons 
absorb  only  a  tiny  fraction  of  the  laser  energy,  about  0.01%.  The  ‘critical  target  thickness’, 
separating  the  two  regimes,  can  be  compared  with  the  skin  depth  Lskin-  At  these  conditions  the 
relativistic  skin  depth  is  "^0.15  /xm,  approximately  one-third  of  the  ‘critical  target  thickness’. 
If  the  target  thickness  greatly  exceeds  the  skin  depth,  i.e.  Lau  ^  Tskin,  the  laser  field  is 
exponentially  attenuated  inside  the  front  layer  and  cannot  penetrate  deep  enough  to  reach  the 
rear  layer,  leading  to  a  decline  of  the  mean  deuteron  energy  and  energy  conversion  efficiency. 
The  position  of  the  fall-off  depends  upon  the  skin  depth,  which,  in  turn,  depends  upon  the 
electron  density.  If  the  electron  density  was  lower,  say  =  lO^c^  the  skin  depth  would  be 
larger  and  the  fall-off  would  occur  for  larger  target  thicknesses.  Our  simulation  results  agree 
well  with  published  data  [41],  which  show  a  rapid  fall-off  of  the  mean  ion  energy  with  the 
target  thickness  increasing.  The  results  in  figures  5(e)  and  (/)  do  not  indicate  an  optimum 
target  thickness  (none  of  the  studied  parameters  goes  through  a  maximum).  The  reason  is 
again  the  assumed  plasma  density  of  the  front  layer.  Studies  have  shown  that  the  optimum 
target  thickness  scales  inversely  with  the  target  plasma  density  [48]  and  in  our  case  the  plasma 
is  highly  overdense  (Hq  =  lOO^c)-  According  to  figure  20(b)  from  [48]  at  ric/riQ  =  0.01 
the  optimum  target  thickness  is  expected  to  be  0.05  /xm,  which  is  below  the  thickness  range 
plotted  in  figure  5.  In  order  to  compare  with  the  data  in  [48]  we  repeated  the  simulations 
with  Hq  =  20nc.  We  observed  a  broad  peak  for  both  the  mean  deuteron  energy  and  energy 
conversion  efficiency  located  between  Lau  =  0.1  and  0.25  /xm,  which  agrees  well  with  the 
value  Lau  =  0.25  /xm  derived  from  the  same  figure.  The  thickness  of  the  deuterium  layer 
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may  also  be  significant.  For  a  given  laser  energy  F^iaser  =  and  energy  conversion 

efficiency  the  mean  ion  energy  is  (wd+)  =  ^D^iaser/A^D,  where  A^d  =  tt/^qLd^d  is  the 
number  of  deuterium  ions.  After  simple  manipulations  one  readily  derives 

{u^+)  =  r]j) — (2) 

The  nominator  is  the  laser  fluence  Fq  =  £iaser/(^^o)  ~  (laser  energy  per  unit  area) 

in  units  J  m“^,  while  the  denominator  contains  the  thickness  of  the  deuterium  layer.  Formula 
(2)  implies  that  the  average  deuteron  energy  depends  on  two  parameters,  the  laser  fluence  and 
the  deuterium  layer  thickness,  provided  the  conversion  efficiency  is  reasonably  constant 
over  a  range  of  parameters.  For  fixed  laser  parameters  this  is  true  for  relatively  thin  targets, 
Lau  ^0.5  /xm,  comparable  to  the  plasma  skin  depth.  Higher  deuteron  energy  can  be  achieved 
by  using  either  a  laser  with  higher  fluence  (decreasing,  if  possible,  the  laser  spot  size  Rq)  or 
a  thinner  deuterium-containing  layer.  The  latter  expresses  the  fact  that  the  energy  absorbed 
by  the  deuterium  layer  is  spread  over  fewer  deuterons,  allowing  them  to  gain  higher  mean 
average  energy.  Indeed,  simulations  with  thickness  of  the  deuterium  layer  Ld  =  0.1,  0.05, 
0.02  and  0.01  /xm  confirm  the  upward  trend  of  the  mean  deuteron  energy  with  Ld  decreasing 
(figure  5(f)).  It  should  be  kept  in  mind,  however,  that  there  is  a  drawback  of  making  the 
deuterium  layer  too  thin;  it  results  in  acceleration  of  fewer  deuterons. 

The  third  and  fourth  columns  in  figure  5  plot  the  same  parameters  versus  peak  laser 
intensity  and  laser  pulse  duration,  respectively.  As  expected,  both  the  energy  conversion 
efficiency  and  the  mean  energy  of  increase  with  peak  laser  intensity  and  pulse  duration. 
For  peak  laser  intensities  approaching  Iq  =  10^^  W  m~^  the  conversion  efficiency  saturates 
and  the  mean  deuteron  energy  increases  primarily  due  to  the  increase  of  the  peak  laser  intensity. 
The  conversion  efficiency  versus  laser  pulse  durations  shows  two  distinct  regimes  (figure  5(/)). 
It  increases  exponentially  for  tq  <  100  fs,  indicating  that  the  time  is  not  sufficient  to  accelerate 
the  ions.  For  tq  >  lOOfs  the  moderate  increase  with  tq  suggests  that  the  ions  have  been 
sufficiently  accelerated  to  their  maximum  velocity.  Very  long  laser  pulse  duration,  however,  is 
to  be  avoided  in  order  to  prevent  Au  ions  from  being  accelerated,  aiming  at  depositing  a  larger 
fraction  of  the  laser  energy  into  the  light  ions. 

Figures  6  to  9  address  the  main  deuteron  beam  parameters,  the  energy  and  angular 
distribution  functions  of  D"^.  They  were  studied  by  varying  the  same  laser  and  target  parameters 
(one  at  a  time)  as  in  figure  5.  The  effect  of  the  target  thickness  is  seen  in  figure  6.  The  energy 
distribution  function  is  bi-Max wellian,  having  a  bulk  of  cold  ions  and  a  long  tail,  which  extends 
to  kinetic  energy  of  several  MeV.  For  Lau  =  0.1  /xm  ions  with  kinetic  energy  of  up  to  15  MeV 
are  observed.  With  increasing  target  thickness  the  maximum  ion  energy  decreases.  The 
decrease  is  most  dramatic  for  target  thickness  between  1  and  3  /xm.  The  angular  distribution 
function  of  is  also  target  thickness  sensitive.  For  Lau  <  1  /^m  the  angular  distribution 
function  peaks  at  very  small  scattering  angle.  With  increasing  target  thickness  the  scattering 
becomes  more  isotropic.  From  figures  5  and  6  we  conclude  that  very  thin  targets,  of  the  order  of 
or  less  than  1  /xm  (for  Hq  =  lOOric),  favour  the  generation  of  energetic  and  forward  accelerated 
ions.  The  thickness  of  the  deuterium  layer  has  little  impact  on  the  deuteron  distribution 
function  (figure  7,  left).  It  has  a  profound  effect,  however,  on  the  angular  scattering  of  the 
deuterons.  The  thinner  the  deuterium  layer,  the  narrower  the  angular  distribution  becomes 
(figure  7,  right).  Therefore,  it  is  advantageous  to  prepare  targets  with  a  very  thin  (10-50  nm) 
layer  of  deuterium-containing  material. 

The  next  sequence  of  figures  illustrates  the  ion  acceleration  versus  input  laser  parameters. 
The  deuteron  energy  and  angular  distribution  function  for  four  peak  laser  intensities  varying 
from  /o  =  10^^  to  /q  =  10^^  W  m~^  are  plotted  in  figure  8.  As  expected,  the  maximum 
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Figure  6.  Deuteron  energy  distribution  function  (left  column)  and  angular  distribution  function 
(right  column)  for  different  target  thicknesses.  All  parameters,  except  the  target  thickness,  are  as 
in  figure  2. 
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ion  energy  increases  with  peak  laser  intensity,  following  approximately  a  linear  dependence. 
For  peak  laser  intensity  equal  to  or  below  Iq  =  10^^  W  m~^  the  ion  distribution  function  is 
bi-Max wellian.  But  for  the  highest  peak  laser  intensity  of  Iq  =  10^^  W  m~^  an  ion  beam  with 
broad  energy  spectrum  between  2  and  50  MeV  is  formed.  The  laser  pulse  duration  affects  only 
the  ion  energy,  but  not  the  angular  distribution  (figure  9).  The  ions  attain  higher  energy  with  the 
pulse  duration  increasing  simply  because  more  laser  energy  is  deposited  and  because  the  ions 
are  accelerated  for  a  longer  period  of  time.  These  results  naturally  pose  the  question,  whether 
the  deuteron  acceleration  depends  on  the  peak  laser  intensity  and  pulse  duration  separately,  or 
on  the  combination  of  both.  We  chose  as  universal  parameter  the  laser  fiuence  Fq  =  IqToII, 
following  the  derivation  of  equation  (2).  In  figure  10  we  plot  the  deuteron  mean  energy  and 
conversion  efficiency  versus  laser  fiuence.  The  solid  curve  refers  to  variation  of  the  peak  laser 
intensity,  while  the  dashed  curve  with  symbols  refers  to  variation  of  the  laser  pulse  duration. 
Not  only  do  both  curves  exhibit  similar  trends  but  they  are  also  quantitatively  similar.  We  can 
conclude  that  the  properties  of  the  deuteron  beam  depend  on  the  laser  fiuence  rather  than  the 
input  laser  intensity  and  pulse  duration  separately. 
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Figure  7.  Deuteron  energy  distribution  function  (left  column)  and  angular  distribution  function 
(right  column)  for  different  thicknesses  of  the  D  layer.  All  parameters,  except  the  thickness  of  the 
D  layer,  are  as  in  figure  2. 


3.4.  Comparison  with  other  simulations 

Benchmarks  of  simulation  results  are  always  essential  to  validate  a  model.  A  one-to-one 
comparison  is  highly  desirable,  but  often  difficult  to  do,  especially  at  the  extreme  conditions 
prevalent  here.  Either  the  laser  or  target  parameters  may  not  be  exactly  matched.  In  this  case, 
we  deal  with  deuterium  and  the  data  for  deuteron  acceleration  is  scarce.  The  available  data  are 
mostly  for  uniform  (not  layered)  targets  [33-35].  Therefore,  due  to  the  lack  of  adequate  data, 
we  compare  deuteron  with  proton  acceleration.  This  is  not  too  restrictive  since  the  mechanism 
of  light  ion  acceleration  from  the  rear  surface  of  the  target  via  TNSA  is  very  general  and  can  be 
applied  to  ions  with  comparable  charge- to-mass  ratio.  Proton  acceleration  has  been  extensively 
studied  and  there  is  a  wealth  of  data  published  in  the  literature.  For  example,  Esirkepov 
et  al  [52]  have  done  extensive  parametric  study  of  proton  acceleration  from  a  double-layer 
target  for  laser  and  target  parameters,  which  are  similar  to  our  calculations.  The  thickness 
of  the  front  and  rear  layers  are  the  same,  1  /xm  and  0.06  /xm,  respectively,  and  the  plasma  is 
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Figure  8.  Deuteron  energy  distribution  function  (left  column)  and  angular  distribution  function 
(right  column)  for  different  peak  laser  intensities.  All  parameters,  except  the  peak  laser  intensity, 
are  as  in  figure  2. 

equally  overdense,  =  lOO^c-  The  data  in  figure  3  from  [52]  plot  the  maximum  ion  energy 
versus  input  laser  energy  and  cover  the  variation  of  a  wide  range  of  peak  laser  intensities, 
pulse  duration  and  laser  spot  size.  We  selected  a  subset  of  these  data  for  laser  beam  diameters 
10  /xm  and  25  /xm  (Rq  =  5  and  12.5  /xm).  We  then  converted  the  laser  energy  to  fiuence  by 
dividing  it  to  the  area  of  the  corresponding  laser  spot  (tv  R^).  Our  simulations  (solid  line)  and 
the  data  from  [52]  (symbols)  are  shown  in  figure  11,  where  the  maximum  ion  energy  is  plotted 
versus  laser  fiuence.  The  agreement  is  very  good,  keeping  in  mind  that  our  simulations  are 
for  deuterons,  while  that  in  [52]  are  for  protons.  The  data  by  Esirkepov  et  al  [52]  for  different 
peak  laser  intensities,  pulse  duration  and  laser  spot  size  also  confirm  what  we  found  in  the 
previous  sub-section,  namely  that  for  sub-micrometer  targets  the  main  factor  governing  the 
ion  acceleration  is  the  laser  fiuence. 


4.  Conclusion 

Numerical  simulations  using  a  2D  electromagnetic  PIC  model  of  deuteron  acceleration  in 
high-intensity  laser  thin-foil  interaction  have  been  presented.  Our  study  investigates  the  effect 
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Figure  9.  Deuteron  energy  distribution  function  (left  column)  and  angular  distribution  function 
(right  column)  for  different  laser  pulse  durations.  All  parameters,  except  the  laser  pulse  duration, 
are  as  in  figure  2. 


of  target  thickness,  peak  laser  intensity  and  laser  pulse  duration  on  the  ion  energy  and  angular 
distribution  functions  and  conversion  efficiency  of  laser  energy  into  ion  kinetic  energy  for 
highly  overdense  plasma  =  100.  Our  findings  can  be  summarized  as  follows: 

(i)  With  increasing  peak  laser  intensity  or  laser  pulse  duration  the  average  deuteron  energy 
and  conversion  efficiency  of  laser  energy  into  kinetic  energy  of  the  deuterons  increase. 
Both  the  average  deuteron  energy  and  the  conversion  efficiency  of  laser  energy  into  kinetic 
energy  of  the  deuterons  depend  essentially  on  the  laser  fiuence  rather  than  the  peak  intensity 
and  pulse  duration  separately. 

(ii)  The  average  deuteron  energy  and  conversion  efficiency  of  laser  energy  into  kinetic  energy 
of  the  deuterons  decrease  as  the  target  thickness  increases.  With  the  increasing  target 
thickness  the  angular  scattering  of  the  deuterons  becomes  more  uniform.  Decreasing  the 
target  thickness  improves  the  angular  collimation  of  D"^. 

(iii)  The  thickness  of  the  deuterium  layer  affects  the  angular  spread  of  deuterons.  The  ion 
beam  collimation  improves  with  the  decreasing  deuterium  layer  thickness. 
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Figure  10.  Deuteron  mean  energy  (a)  and  energy  conversion  efficiency  versus  laser  fluence 
(b).  Solid  line — variation  of  the  peak  laser  intensity  for  constant  laser  pulse  duration,  symbols — 
variation  of  the  laser  pulse  duration  for  constant  peak  laser  intensity. 


Figure  11.  Deuteron  and  proton  maximum  kinetic  energy  versus  laser  fluence.  Solid  line — model: 
Xq  =  1 /xm,  Lau  =  1 /xm,  Ld  =  0.06 /xm,  Hq  =  lOO^c,  Rq  =  3/xm.  Symbols:  Xq  =  1 /xm, 
La\x  =  1  =  0.06 /xm,  Hq  =  lOO^c,  Rq  =  5  /xm  (closed  symbols),  Rq  =  12.5  /xm  (open 

symbols)  [52].  The  solid  line  refers  to  D"^,  while  the  symbols  are  for  H"^. 


(iv)  The  deuterons  are  accelerated  in  the  forward  direction  by  TNSA  for  /q  >  10^3  Wm-2 
and  Lau  <  1  /^m. 

The  simulation  results  demonstrate  that  a  well-collimated  relativistic  beam  of  deuterons  can 
be  produced  in  laser-thin  foil  interaction,  which  can  serve  as  a  precursor  for  generation  of  a 
high-flux  neutron  beam.  These  data  will  be  used  as  input  parameters  of  a  model  simulating 
the  production  and  scattering  of  the  neutrons. 
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Appendix 


The  problem  modeled  is  studied  in  a  2D  Cartesian  coordinate  system.  There  is  no  variation 
in  the  z  direction,  i.e.  d/dz  =  0  for  all  variables.  We  consider  a  TE  electromagnetic 
wave  propagating  in  the  +v  direction  with  field  components  E  =  {Ex{x,  y),  Ey{x,  y),0), 
H  =  (0,  0,  Hz{x,  y))  and  B  =  With  these  assumptions  the  Maxwell  equations  takes 

the  forms 


with  j  =  {jx  ,  y ) ,  jy  ,  y ) ,  0)  being  the  conduction  current  density.  The  relativistic  equations 
of  motion  of  a  charged  particle  a  in  an  electromagnetic  field  is  described  by  the  following  set 
of  equations: 


d? 


=  q^\E  + 


Pa 

niaYa 


(2a) 


v„  = 


Pa/ mo, 


\/l  +  (Pa/maC)^ 


(2b) 


AXa 


(2c) 


where  pa  =  (PaA^^y)^  Pa,y(x.y),Q),  V„  =  (Vo,,;t(x,  );),  );),  0),  Ma  and  qoi  are 

the  particle  relativistic  momentum,  velocity,  mass  and  charge,  respectively,  and  = 
y/\  +  {palniac)'^  is  the  relativistic  factor.  Equation  ila)  is  discretized  to  form  two  ED  equations 
for  the  impulse  components  with  a  solution 


n+l 

i^a,x 


qoiET^^^At  +  riT^''^qo,E’Z^'^ At  +  (l  -  (rh^'^y')pl^: 


\  +  (rir'^)^ 


(3a) 


n+l 


q„Ef'^At  -  qTI\„ETI^At  +  (l  -  (pT'^f)pA  -  2qTI^  i 

1  + 


(2>h) 


The  dimensionless  parameter  /2m is  the  cyclotron  frequency 

times  the  time  step.  The  superscript  denotes  quantities  at  the  corresponding  time  level. 
The  current  density  at  the  new  time  level  =  Uo^qo/v/^^  is  calculated  by  multiplying 
(3)  by  the  factor  riaqa/ .  Second-order  correct  discretization  of  equations  (la)-(lc) 
mandates  the  current  density  in  the  Maxwell  equations  to  be  calculated  at  half  time  steps, 
which  we  do  by  averaging  over  4  and  4+1,  expression  for  the  v 

and  y  components  of  ja  takes  the  forms 


•n+1/2  _  ^n+l/2  rn+1/2  Ln+1/2  j^n+l/l  a  -n 

•w+1/2  _  ^«+l/2  rn+1/2  _  i^n+l/2  pn+l/2  a  •« 
Ja,y  ~  ^y  ^x  ^  ^Ja,y^ 


(4a) 

(4b) 
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where 


Ja,x  + 


A/" 

Ja,y 


Hy  -  Ja,x 

1 + 


n+\/2  ^  n^qlAt 

and  .  The  global  current  density  j  =  jc,  is  computed  by  adding  the 

current  densities  of  individual  particles: 


(r )  =  a  (r )  (r )  +  b(r)Ey  (r )  +  Aj^  (r ) ,  (5a) 

jy  (r)  =  a (r) Ey  (r)  -  b (r) E^  (r)  +  Ajy  (r ) .  (5b) 


The  coefficients  are  calculated  on  grid  nodes  r  by  summing  the  coefficients  of  individual 
particles  times  the  particle  shape  function.  The  key  issue  is  to  collect  (and  interpolate)  the 
coefficients,  rather  than  the  current  density  itself.  Equations  (5a)  and  (5b)  are  written  in  such 
a  fashion,  that  the  electric  field  is  factored  out.  Inserting  (5a)  and  (5b)  into  the  Maxwell 
equations  (la)-(lc)  and  transferring  the  terms  with  the  electric  field  into  the  left-hand  side 
we  get 


dE, 

dt 


a  b 

4"  -  Ex  4"  -  Ey 

£o  ^0 


dEy  a  b 

^  +-£v - Ey, 


dt 


^0 
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^0 

1 

^0 


(JJh 

V  dx 


(6a) 

(6b) 


dH, 

^±1 

(dE, 

dEy 

dt 

/^o 

\  dy 

dx 

(6c) 


Solving  (6a)-(6c)  is  advantageous  compared  with  solving  the  original  equations  (la)-(lc). 
The  leading  coefficient  a  is  non-negative  and  changes  on  a  time  scale,  much  longer  than  a 
laser  cycle.  The  system  of  two  linear  equations  for  Ex  and  Ey  has  a  determinant  that  is  always 
positive  due  to  the  particular  arrangement  of  the  coefficients  a  and  Z?.  As  a  consequence,  the 
system  of  equations  (6a)-(6c)  is  very  robust  and  efficient  in  dumping  spurious  solutions. 

Equations  (6a)-(6c)  require  initial  and  boundary  conditions.  At  the  beginning  of  the 
computations  the  electromagnetic  field  components  are  set  to  zero.  To  allow  outgoing 
electromagnetic  waves  to  leave  the  computational  domain  without  spurious  reflections  from 
the  edges,  we  use  boundary  conditions  [69] 


d  d  " 

dt  dx  ^ 

jr  =  0, 

V  =  0, 

{Id) 

d  d  \ 

dt'^^^J 

1  T  =  0, 

—  Tj., 

{lb) 

1 

1 

1  y  =  0, 

y  = 

(7c) 

d  d  \ 

dt  dyj 

1  y  =  0, 

y  —  ^y/2, 

{Id) 

where  Y  =  Ex,  Ey.  There  are  no  boundary  conditions  for  :  it  is  calculated  directly  from 
(6c).  Equations  (6a)-(6c)  also  require  a  ‘source  field’  located  at  the  left  boundary  v  =  0.  Eor 
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each  time  step  boundary  conditions  ila)-{ld)  are  applied  and  then  Ey  is  updated  according 
to  [70]: 


Ey{y,  t)  Ey(y,  t)  —  IcoQ/ytE^iy,  t)  sin(a;oO.  (8) 

where  E{){y,t)  =  y/2I(y,  t)/c£o  is  the  envelope  of  the  unperturbed  laser  electric  field. 
The  source  electric  field  is,  in  fact,  twice  the  time  derivative  of  the  driver  laser  field 
Eo(y,  t)  cos(a;oO-  Note  that  in  this  discretization  scheme  the  field  components  and  the  current 
density  are  located  on  grid  nodes,  which  is  different  from  the  arrangement  in  the  conventional 
Yee  algorithm. 
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Abstract 

The  angular  distribution  of  neutrons  formed  in  nuclear  fusion  reactions  of 
a  high-energy  deuteron  beam  with  a  deuterated  polyethylene  (CD2)  was 
investigated  with  a  Monte  Carlo  ion  beam-target  deposition  model.  The  initial 
conditions  were  obtained  from  a  two-dimensional  particle-in-cell  laser-target 
deposition  model.  The  neutron  yield  and  its  angular  distribution  were  studied 
as  a  function  of  peak  laser  intensity,  laser  pulse  duration  and  primary  target 
thickness.  The  proposed  scheme  for  neutron  production  delivers  a  typical 
neutron  yield  of  10“^-10“^  neutrons/ion  and  10^-10^  neutrons  laser  energy. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

High-intensity  ultrashort  pulse  lasers  are  compact  and  versatile  systems  that  can  be  employed 
to  induce  nuclear  fusion  reactions.  If  laser  radiation  of  sufficiently  high  intensity  irradiates 
a  thin  planar  foil,  fast  ions  with  energies  in  the  megaelectronvolt  range  are  generated  on 
the  back  surface  of  the  foil.  These  ions  can  drive  nuclear  fusion  reactions  and  generate 
a  copious  amount  of  neutrons.  Such  laser-driven  neutron  sources  are  an  alternative  to  the 
accelerator-  and  reactor-driven  sources  offering  high  brightness,  compactness,  short  duration 
and  relatively  low  cost.  Applications  include  fast  neutron  radiography  [1],  transmutation  of 
nuclear  waste  [2]  and  fusion  research  [3,4].  Another  exciting  application,  which  is  currently 
under  intense  investigation,  is  the  generation  of  high  brightness  neutron  beams  for  neutron 
activation.  Figure  1(a)  is  a  schematic  visualization  of  the  neutron  activation  process.  During 
the  first  stage,  a  collimated  high-energy  ion  beam  (~MeV)  is  generated  from  the  primary 
target.  In  the  next  stage,  fast  neutrons  are  born  in  nuclear  fusion  reactions  driven  by  the 
fast  ions  from  a  thick  (^  1  mm)  secondary  target  made  of  deuterated  material.  During  the 
third  stage,  the  neutrons  impinge  on  the  target  of  interest  and  the  interaction  of  the  fast 
(several  megaelectronvolt)  neutrons  with  the  target  leads  to  the  emission  of  gamma  rays.  The 
gamma-ray  ‘signature’  is  unique  for  each  element  and,  in  principle,  can  be  used  to  identify 
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Figure  1.  Neutron  production  concept  (a)  and  neutron  activation  of  a  distant  target  (b).  The  laser 
pulse  propagates  from  left  to  right  along  the  axis  ‘x’. 


the  constituents  of  the  material  producing  them.  The  most  challenging  aspect  of  the  neutron 
activation  process  is  to  deliver  a  sufficient  number  of  neutrons  on  targets  located  ~5-50  m 
from  the  neutron  source  (figure  1(b)).  The  problem  is  that  neutrons  born  in  nuclear  fusion 
reactions  are  scattered  (more  or  less)  isotropically  and  only  a  tiny  fraction  of  the  total  number 
of  neutrons  will  hit  a  distant  target.  Thus  the  problem  is  not  only  to  produce  a  large  number 
of  neutrons  but  also  to  scatter  them  preferentially  in  the  forward  direction.  In  fact,  producing 
smaller  numbers  of  well  collimated  neutrons  may  be  advantageous  compared  with  producing  a 
large  number  of  isotropically  scattered  neutrons  for  two  reasons.  First,  the  number  of  neutrons 
on  target  is  determined  by  the  neutron  flux  and  not  by  the  number  of  neutrons,  and  second, 
shielding  of  personnel  and  equipment  may  be  alleviated  or  eliminated.  One  of  the  options  is 
to  use  relatively  inexpensive  table-top  high-intensity  lasers,  which  potentially  can  generate  a 
collimated  beam  of  neutrons. 

Neutron  production  from  high-intensity  laser-target  interaction  has  already  been 
demonstrated  in  a  number  of  experimental  [1, 5-17]  and  theoretical  [3, 14-22]  studies.  The 
bulk  of  the  work  has  been  done  during  the  last  ten  years,  when  ultrashort  pulse  lasers 
with  parameters  suitable  for  neutron  production  became  available.  Most  experiments  were 
performed  with  lasers  having  energies  ranging  from  a  few  Joules  to  about  100  J.  The  total 
neutron  yield  reported  is  typically  10"^- 10^  neutrons.  Due  to  the  differences  in  laser  systems,  in 
particular,  delivered  laser  energy  on  target,  it  is  better  to  compare  the  neutron  yield  per  unit  laser 
energy.  Still,  the  results  vary  widely:  from  ~  10"^  neutrons  J“H7, 8, 12]  to  ~  10^  neutrons  J^H^] . 
A  summary  of  the  neutron  yield  from  various  experiments  is  given  in  table  1 .  The  angular 
distribution  of  neutrons  is  of  critical  importance  for  producing  a  high  flux  of  neutrons.  Toupin 
et  al  studied  the  angular  distribution  of  neutron  emission  as  a  function  of  laser  intensity  and 
maximum  electron  density  [17].  Their  results  indicate  that  if  the  medium  generating  light  ions 
is  highly  overdense  plasma  (such  as  from  laser-thin-foil  interaction),  the  neutrons  are  emitted 
preferentially  in  the  forward  direction.  The  neutron  production  model  employed  was  based 
on  the  continuous  slowing  down  approximation  of  deuterons  in  a  solid  target.  The  deuteron 
constantly  loses  energy  until  its  energy  drops  below  a  prescribed  cut-off  value.  The  neutron 
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Table  1.  Laser  and  target  parameters  and  neutron  yield  from  high-intensity  laser  solid  interactions. 


Target 

£(J) 

Iq  (Wcm-2) 

Tq  (ps) 

Y 

r(nJ-‘) 

References 

CgDg 

20 

8  X  10‘* 

1.3 

1  X  10’  (nsr-‘) 

4  X  10’ 

[5] 

CD2 

2 

3  X  10‘® 

0.3 

7  X  10'’  (n) 

4  X  10'’ 

[6] 

CD2 

20 

3  X  10‘* 

0.45 

9  X  10^  (n) 

5  X  10'* 

[7] 

CD2 

0.3 

2  X  10‘* 

0.05 

1  X  lO'*  (n) 

3  X  10'* 

[8] 

C2D4 

0.2 

1  X  10‘* 

0.16 

1.4  X  10^  (n) 

7  X  10^ 

[10] 

CD 

90 

5  X  10‘® 

0.9-1. 2 

3  X  10'’  (nsr-') 

4  X  10^ 

[11] 

CD2 

5 

2  X  10‘* 

1.5 

1  X  10^  (n) 

2  X  10'* 

[12] 

CD 

20-50 

2  X  10‘® 

0.45-0.8 

1  X  lO’'  (n) 

5  X  10^ 

[15] 

CgDg 

2.3 

2  X  10‘® 

0.5 

7  X  10"*  (nsr“*) 

2  X  10'* 

[16] 

yield  is  calculated  via  integration  of  the  neutron  production  differential  cross  section  over  the 
deuteron  energy  and  angular  distribution  functions.  One  of  the  main  simplifications  of  the 
model  is  that  the  deuteron  trajectory  is  a  straight  line  (but  as  we  will  see  in  section  4,  this  is 
actually  a  very  good  approximation).  To  overcome  this  limitation,  Habara  etal  [15]  and  Izumi 
et  a/  [16]  used  a  3D  Monte  Carlo  code,  which  tracks  the  ion  trajectory  through  the  secondary 
target. 

Our  previous  study  focused  on  the  first  stage  of  neutron  production,  i.e.  the  formation  of  a 
high-energy  deuteron  beam  from  high-intensity  laser-target  interactions  [23] .  The  laser-target 
deposition  model  of  the  first  stage  is  a  two-dimensional  electromagnetic  particle-in-cell  model. 
It  was  employed  to  study  theoretically  the  production  of  megaelectronvolt  deuterons  during  the 
interaction  of  intense  ultrashort  laser  pulses  with  a  double-layer  thin  foil.  The  directionality 
and  energy  of  the  deuteron  beam  were  studied  as  a  function  of  peak  laser  intensity,  laser  pulse 
duration  and  target  thickness.  This  paper  focuses  on  the  second  stage  of  neutron  generation. 
An  ion-beam-target  deposition  model  is  developed  to  calculate  the  neutron  production.  The 
ion  transport  through  a  thick  deuterated  target  (CD2)  is  simulated  using  three-dimensional 
Monte  Carlo  calculations,  as  in  [15, 16].  Each  simulation  yields  the  angular  distribution  of 
neutrons  and  the  total  neutron  yield  in  absolute  units  as  a  result  of  the  nuclear  fusion  reaction 
D-i-D  ^  ^He-i-n.  The  model  is  used  to  perform  a  comprehensive  parametric  theoretical  study 
of  the  neutron  generation,  which  includes  a  variation  of  the  peak  laser  intensity,  laser  pulse 
duration  and  thickness  of  both  layers  of  the  primary  target.  The  neutron  yield  in  absolute  units 
and  its  angular  distribution  are  presented  as  a  function  of  these  parameters.  Section  2  details 
the  ion  beam-target  deposition  model.  In  section  3  we  present  the  primary  and  secondary 
target  configurations.  A  simple  model  for  neutron  production  with  forward  directed  mono- 
energetic  ions  and  a  more  detailed  model,  which  incorporates  realistic  deuteron  energy  and 
angular  distributions,  are  compared  in  section  4.  In  this  section  we  also  perform  a  parametric 
study  of  the  angular  flux  of  neutrons,  and  summarize  our  findings  in  section  5. 

2.  Ion  beam-target  deposition  model 

2.7.  Monte  Carlo  simulations 

We  start  this  section  by  briefly  describing  the  output  of  the  laser-target  deposition  model  in 
order  to  explain  how  it  transitions  to  the  ion  beam-target  deposition  model  described  below. 
For  given  laser  and  primary  target  parameters  at  the  end  of  the  simulation  the  deuteron  energy 
and  angle  with  respect  to  the  target  normal  (the  -Fx  direction)  are  recorded  and  stored  in  a  file. 
Both  parameters  are  known  from  the  particle  movement  module  of  the  electromagnetic  model. 
These  data  serve  as  input  of  the  ion  beam-target  deposition  model. 
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The  ion  beam-target  deposition  model  closely  follows  the  neutron  production  models 
developed  by  Habara  et  al  [15],  Izumi  et  al  [16]  and  Toupin  et  al  [17].  In  essence,  this  is  a 
3D  Monte  Carlo  code,  which  follows  the  transport  of  ions  through  the  secondary  target.  Ions 
are  launched  one  by  one  from  the  surface  of  the  secondary  target  with  a  given  energy  and 
angle  01^^  with  respect  to  the  normal  surface.  The  information  for  each  ion  comes  from  the 
output  file  of  the  laser-target  deposition  model,  as  described  in  the  previous  paragraph.  The 
ion  transport  through  the  secondary  target  is  based  upon  the  following  assumptions: 


•  The  secondary  target  is  sufficiently  thick  (a  few  millimeters)  so  that  even  the  most  energetic 
deuterons  are  completely  stopped  inside.  Its  normal  points  in  +x,  the  direction  of  laser 
propagation. 

•  The  ion  distribution  is  axisymmetric  with  respect  to  the  laser  propagation  direction. 

•  The  neutron  distribution  is  also  axisymmetric  with  respect  to  the  same  axis. 

•  The  ion  transport  is  based  on  the  continuous  slowing  down  approximation  using  the 
concept  of  stopping  power  and  angular  scattering  via  small-angle  scattering  in  collisions 
with  the  atoms  of  the  secondary  target.  The  stopping  power  depends  on  the  incident  ion 
energy,  target  material  and  target  density. 

The  ion  transport  through  the  secondary  target  is  modeled  by  following  the  ion  energy 
and  trajectory  in  space  according  to  the  following  algorithm.  The  ion  is  assigned  initial  spatial 
coordinates  g  =  (0,  0,  0),  energy  Ei  =  Ef^^  and  a  velocity  vector  defined  by  a  polar  angle 
0-^  =  with  respect  to  the  normal  surface  and  an  azimuthal  angle  cpi  chosen  randomly  in  the 
interval  [0-27r].  The  ion  advances  in  steps  traversing  a  distance  A/,  which  is  chosen  according 
to  the  initial  ion  energy  so  that  the  deuteron  makes  a  sufficient  number  of  steps  ("^10^  or  more) 
until  it  stops.  The  magnitude  of  the  energy  loss  per  step 

A^i  =  S{Ei)/^l  (1) 


is  determined  by  the  stopping  power  S{Ei).  The  deuteron  is  then  moved  to  a  new  position 
G  ^  G+ A?!,  Ari  =  (A/  sin  6i  cos  ^i.  A/  sin  Oi  sin  cpi,  Al  cos  ^i),  where  it  experiences  apseudo- 
scattering  and  changes  direction  by  selecting  new  azimuthal  and  polar  angles.  The  azimuthal 
angle  cp^^  in  the  center-of-mass  system  is  chosen  randomly  in  the  interval  [0-27r].  The  polar 
angle  0^^  in  the  same  system  is  sampled  from  the  multiple- scattering  distribution  function  for 
the  projected  angle  of  scattering 

(-!))■  ® 

which  is  appropriate  for  describing  multiple  small-angle  scattering  of  high-energy  ions  in 
collisions  with  the  target  atoms  of  the  secondary  target.  A  deuteron  traversing  a  distance 
Al  is  scattered  at  an  angle  whose  mean  square  {0^)  is  given  by  the  well-known  Jackson’s 
formula  [24] 


{e^}  =  47tNt 


/  4Zte^ 
V4:7reo£’i 


2 

ln(204Zj"‘^^)A/, 


(3) 


where  Nt  and  Zt  are  the  secondary  target  density  and  charge  number,  respectively  (e  is  the 
electron  charge  and  Sq  is  the  permittivity  of  free  space).  The  conversion  of  the  post-collision 
direction  of  the  deuteron  ^2^”^  =  <pf^}  from  the  center-of-mass  system  to  the  laboratory 

system  ^2^  =  {^i,  (pi)  determines  the  new  deuteron  direction  in  the  laboratory  system  after  the 
pseudo-collision  takes  place  [25].  The  above  algorithm  is  applied  at  each  step  until  the  ion 
energy  becomes  less  than  the  cut-off  energy  KeV). 
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Figure  2.  Differential  cross  section  for  neutron  production  da(Fi,  6>n)/df^n  for  deuteron  energies 
Ei  of  0.1  MeV  (a),  0.2  MeV  {b),  0.5  MeV  (c),  1  MeV  (J),  2MeV  {e)  and  5  and  lOMeV  (/).  Data 
from  [26]. 


During  each  step  the  neutron  flux  to  speciflc  directions  of  observation  is  also  calculated. 
The  increment  of  the  neutron  flux  at  polar  angle  into  the  solid  angle  d^2n  for  deuteron 
traversing  a  distance  A/  is  calculated  according  to 


(dYM\ 

V  ) 


da(£i,^n) 

d^2n 


(4) 


where  fxo  is  the  density  of  deuterium  substance  in  the  secondary  target  (a  value  provided  in 
section  3)  and  da(Ei,  ^n)/d^2n  is  the  differential  cross  section  for  neutron  production^  [26], 
plotted  in  figure  2.  It  is  important  to  note  that  the  neutron  production  direction  in  equation  (4) 
is  accounted  for  in  a  coordinate  system  with  axis  pointing  in  the  direction  of  the  deuteron 
velocity,  which  is  rotated  with  respect  to  the  laboratory  coordinate  system.  In  other  words, 
is  measured  with  respect  to  Oi,  while  the  polar  angles  c^n  and  Oi  are  measured  (in  the 
laboratory  system)  with  respect  to  the  secondary  target  normal.  A  transformation  of  the 
directional  vector  from  the  rotated  coordinate  systems  to  that  of  the  laboratory  system  yields 
0^  =  cos“^  (cos  Oi  cos  «n  +  sin  Oi  sin  an  cos  cpi),  which  directly  relates  angle  an  to  ^n  via  the  ion 
directional  angles  {Oucpi}  [17]. 


2.2.  Simple  model:  mono -energetic  deuteron  beam  normally  incident  on  a  target 

It  is  fairly  straightforward  to  calculate  the  neutron  yield  in  the  special  case  when  the  deuterons 
are  mono-energetic,  normally  incident  upon  the  secondary  target  and  move  in  a  straight  line 
inside  the  secondary  target  until  they  stop.  The  neutron  yield  at  an  angle  0^  =  ^  0^  ^  7t) 

relative  to  the  direction  of  the  deuteron  beam  is  calculated  by  integrating  the  ratio  of  the  nuclear 

^  The  cross  section  is  tabulated  only  for  incident  deuteron  energies  below  13.8  MeV.  For  Ei  >  13.8  MeV  we  scaled 
the  magnitude  of  the  total  cross  section  and  assumed  the  same  angular  distribution  as  that  for  the  largest  tabulated 
energy  available  (13.8  MeV). 
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fusion  differential  cross  section  to  the  stopping  power  from  the  initial  deuteron  energy  to 
the  final  cut-off  energy  £;pt-off 


(5) 


Formula  (5)  is  the  number  of  neutrons  per  steradian  per  incident  ion.  In  spite  of  the 
simplifications  involved,  formula  (5)  is  very  useful.  The  neutron  yield  can  be  tabulated  for  a 
set  of  incident  ion  energies  and  neutron  direction  angles  6^.  Further  integration  over  the 
ion  energy  distribution  yields  the  total  neutron  production,  as  was  done  in  [17]. 

3.  Primary  and  secondary  target  parameters 

The  primary  target  is  a  two-layer  thin  foil  with  thickness  between  0.1  and  1  /xm.  The  front 
layer  is  made  of  high-Z  material.  The  rear  layer  is  very  thin  (0.01-0.1  /xm)  and  is  made 
of  deuterium-containing  material.  The  choice  of  high-Z-low-Z  material  is  motivated  by  the 
observation  that  if  the  ratio  of  mass  to  charge  for  the  front  and  rear  material  is  sufficiently  large, 
the  light  ions  from  the  rear  surface  are  accelerated  much  more  efficiently  than  the  heavy  ones. 
The  front  layer  material  is  gold,  a  common  substance  in  this  kind  of  high-intensity  laser-target 
interactions.  The  plasma  slab  corresponding  to  this  layer  has  a  thickness  Lau,  density  ^au  and 
charge  Zau  —  1,  respectively.  For  simplicity  we  assumed  that  the  rear  layer  is  made  of  pure 
deuterium  with  thickness  and  density  Ld  and  respectively.  Both  the  front  and  the  rear 
layers  have  the  form  of  a  disc  with  a  prescribed  diameter  Dau  and  Dd,  respectively. 

The  secondary  target  is  a  thick  slab  of  deuterated  material.  We  chose  deuterated 
polyethylene  (CD2),  which  is  commonly  used  in  both  experiments  and  simulations.  The 
density  of  CD2  (1.105  x  10^kgm“^)  [27]  determines  the  fractional  densities  of  carbon  and 
deuterium  in  the  target:  =  4.1  x  10^^  m“^  and  =  8.2  x  10^^  m“^,  respectively.  The 

material  and  its  density,  in  turn,  determine  the  stopping  power  of  deuterons  in  the  secondary 
target.  Since  the  stopping  power  depends  on  the  projectile  velocity,  the  stopping  power 
of  in  carbon  and  deuterium  is  computed  by  shifting  the  corresponding  stopping  power 
for  in  the  energy  domain,  i.e.  =  S’f  “^(£^/2),  where  Z  =  C,  D  stands 

for  carbon  and  deuterium  [28,29].  The  computed  individual  stopping  powers  per  atom 
~^(£^)  and  5'}^  “^(£)  are  shown  in  figure  ?>{a)  and  the  total  stopping  power  of  in  CD2 
5'D^-cd2(^)  =  sY~^{E)n'l^  is  shown  in  figure  The  latter  is  the  stopping 

power  in  the  right-hand  side  of  equation  (1). 

4.  Results  and  discussions 

4.1.  Neutron  yield  from  a  mono -energetic  deuteron  beam  normally  incident  on  a  target 

We  commence  with  the  simplest  case  of  a  mono-energetic  beam  of  deuterons  normally  incident 
upon  the  secondary  target  (section  2.2).  We  performed  a  series  of  simulations  by  launching 
deuterons  with  energies  Ef^^  =  0.1,  0.3,  1,  3,  10  and  30MeV  at  incident  angle  normal  to 
the  surface.  The  angular  distribution  of  the  neutron  yield  per  deuteron  for  various  incident 
deuteron  energies  is  shown  in  figure  4  (solid  lines).  It  can  be  separated  into  two  regions, 
depending  on  the  incident  ion  energy:  for  deuteron  energy  below  about  1  MeV  the  angular 
distribution  is  close  to  uniform,  while  for  higher  energies  it  is  peaked  forward.  The  total 
neutron  yield  (integrated  over  the  scattering  angle),  together  with  the  total  cross  section  for 
neutron  production,  is  plotted  in  figure  5.  The  total  neutron  yield  increases  nonlinearly  with 


6 


Plasma  Phys.  Control.  Fusion  50  (2008)  065016 


J  Davis  and  G  M  Petrov 


E 

> 


in 


O 


0) 

o 

Q. 

O) 

c 

■q. 

Q. 

O 

CO 


Figure  3.  Stopping  power  per  atom  of  in  C  and  D  (a)  and  total  stopping  power  of  in  CD2 
(b)  versus  deuteron  energy.  Symbols:  stopping  power  per  atom  of  D+  in  C  [34]. 


e  (deg) 

Figure  4.  Angular  distribution  of  the  neutron  yield  per  ion  0n)/df^  versus  scattering 

angle  for  various  incident  deuteron  energies.  Simple  model — solid  lines,  symbols — Monte  Carlo 
model. 


the  deuteron  energy  and  can  be  well  approximated  as 

[7.8  X  ^  IMeV, 

y/rinc\  _  I  ^  1  ^  1 

‘  [7.8  X  10-®(£'‘“)2-2  ^  IMeV. 

The  near-quadratic  dependence  of  the  neutron  yield  versus  incident  deuteron  energy  suggests 
that  in  simulations  with  deuterons  having  sufficient  energy  spread  the  high-energy  deuterons 
from  the  tail  of  the  ion  energy  distribution  function  (lEDF)  may  contribute  significantly  to  the 
total  neutron  yield.  This  point  is  proven  in  the  following  sub-section. 

Our  next  task  is  to  study  the  impact  of  scattering  on  the  neutron  yield.  We  performed 
another  series  of  simulations  by  launching  mono-energetic  ions  with  the  same  incident  energy 
and  angle.  The  only  difference  from  the  previous  set  of  simulations  is  that  now  the  deuteron 
trajectory  is  modeled  with  Monte  Carlo,  rather  than  being  a  straight  line.  For  improved 
statistics,  the  neutron  yield  is  averaged  over  ~10^  ions  (for  each  incident  ion  energy).  The 
angular  distribution  of  the  neutron  yield  per  ion  is  plotted  in  figure  4  (symbols).  The  results 
show  that  the  actual  deuteron  trajectory  is  of  little  importance  for  the  neutron  yield.  This  is 
confirmed  in  figure  6,  where  sample  deuteron  trajectories  are  plotted  for  the  same  incident 
deuteron  energies.  The  abscissa  is  the  deuteron  position  in  the  direction  parallel  to  the  angle 
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Figure  5.  Total  cross  section  for  neutron  production  cr(Fi)  =  (dcr(£i,  0n)/d^n)  (a)  and 

total  neutron  yield  Y {E™^)  =  (dF 6>n)/df2n)df^n  versus  deuteron  energy  (b). 
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Figure  6.  Deflection  in  the  transverse  direction  y/y'^  -\-z^  versus  longitudinal  position  of  deuteron 
trajectories  for  incident  deuteron  energies  0.1  MeV  (a),  0.3  MeV  {h),  1  MeV  (c),  3MeV  (J), 
lOMeV  {e)  and  30 MeV  (/). 

of  incidence  (laser  propagation  direction  +x),  while  the  ordinate  plots  the  deflection  -yjy^  + 
in  the  transverse  direction.  The  transverse  deflection  at  the  end  of  the  deuteron  trajectory  is 
very  small  compared  with  the  penetration  depth  and  only  after  the  deuteron  has  lost  most  of 
its  initial  energy.  Since  the  bulk  of  the  neutron  production  occurs  at  high  deuteron  energy,  the 
last  segment  of  the  deuteron  trajectory  has  negligible  contribution  to  the  total  neutron  yield. 
These  results  confirm  that  if  the  deuteron  trajectory  is  assumed  to  be  a  straight  line,  the  neutron 
production  is  accurately  computed  in  this  regime  of  laser-target  interactions. 


4.2.  Parametric  study:  variation  of  the  peak  laser  intensity,  laser  pulse  duration  and  primary 
target  thickness 

In  this  sub-section  a  parametric  study  of  the  angular  distribution  of  neutrons  and  neutron 
yield  (in  absolute  units)  is  performed.  Both  are  calculated  as  a  function  of  peak  laser 
intensity,  laser  pulse  duration  and  target  thickness.  The  laser  and  primary  target  parameters 
are  varied  in  a  range  suitable  for  neutron  generation:  peak  laser  intensity  lO^^-lO^^  W  m“^ 
(lO^^-lO^^  W  cm“^),  deuterium  layer  thickness  0.01-0.1 /xm  and  gold  layer  thickness 
0.1-1  /xm.  The  laser  radiation  is  an  ultrashort  pulse  of  linearly  polarized  light  normally 
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Figure  7.  Deuteron  energy  distribution  function  (a)-(d),  angular  distribution  function  (e)-(h), 
neutron  yield  per  ion  versus  scattering  angle  (/)-(/)  and  total  neutron  yield  versus  scattering 
angle  (m)-(p)  for  different  thicknesses  Lau  of  the  Au  layer.  The  laser  and  target  parameters 
are  Iq  =  10^^Wm“^,  Xq  =  1 /xm,  Rq  =  3 /xm,  to  =  80  fs,  ^laser  =  1.13  J,  Ld  =  0.05 /xm, 
Uq  =  lOO^c. 

incident  on  the  target,  propagating  in  the  +x  direction.  The  laser  pulse  has  a  spatio-temporal 
profile  I (y,t)  =  Iq  sin^(7r  t/zo)  exp(— (y/7?o)^)  with  Iq  being  the  peak  laser  intensity,  Tq  the 
laser  pulse  duration  and  Rq  the  laser  spot  radius.  Details  of  the  simulation  parameters  are 
given  in  [23]. 

The  angular  distribution  of  the  neutron  yield  per  ion  and  the  total  neutron  yield  for  different 
thicknesses  of  the  Au  layer  is  shown  in  figure  7.  The  deuteron  angular  and  energy  distribution 
function,  which  have  been  analyzed  in  [23],  are  shown  for  reference.  The  laser  and  target 
parameters  are  peak  laser  intensity  Iq  =  lO^'^Wm"^,  laser  wavelength  Xq  =  I  /xm,  laser 
spot  size  Rq  =  3  /xm,  laser  pulse  duration  tq  =  80  fs  (tfwhm  =  To/2  =  40  fs),  input  laser 
energy  £^iaser  =  ^^o^qi^fwhm  =  1.13  J  and  deuterium  layer  thickness  Ld  =  0.05  /xm.  The 
primary  target  is  an  overdense  plasma  with  Hq  =  lOOfXc-  The  angular  distribution  of  neutrons 
is  computed  in  elementary  angular  bins  =  (n  l/2)AQ'n,  n  =  0,  . . . ,  35.  The  width  of 
each  angular  bin  is  A^n  =  7r/36  =  8.73  x  10“^  rad  (corresponding  to  5°).  The  neutron  yield 
(the  average  number  of  neutrons  per  ion  per  radian)  to  the  directions  of  observation 

Yr,(an)  =  27t  sm(Q'n) ^77;^ —  (6) 

d^2n 
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is  plotted  in  figures  7  (/)-(/).  The  dashed  line  refers  to  the  simplistic  model  from  the  previous 
sub- section  (mono-energetic  deuteron  beam  with  incident  energy  equal  to  the  mean  deuteron 
energy  and  normally  incident  on  the  target).  The  Monte  Carlo  calculations  differ  from  the 
simplistic  ones  in  two  respects:  the  lEDF  (mono-energetic  versus  actual)  and  the  ion  angular 
distribution  function  (normal  incidence  angle  versus  actual  angle  of  incidence).  The  latter 
plays  a  minor  role  as  the  deuteron  beam  is  only  slightly  divergent,  as  seen  in  the  second  row 
of  figures.  The  major  difference  is  expected  to  be  the  impact  of  the  lEDF.  As  discussed  in 
the  previous  sub-section,  the  neutron  yield  increases  nonlinearly  with  the  deuteron  energy.  In 
addition,  the  lEDF  has  a  very  long  tail  with  maximum  (cut-off)  deuteron  energy  ~10  times 
the  mean  deuteron  energy  (figure  7,  first  row  of  figures).  The  third  row  of  figures  provides  a 
direct  comparison  between  the  simplistic  and  actual  cases.  The  simplified  version  of  neutron 
production  underestimates  the  actual  neutron  yield  by  a  factor  of  two.  Not  surprisingly,  the 
angular  distribution  of  the  neutron  yield  in  both  cases  is  similar,  since  most  of  the  ions  have  a 
very  small  angle  of  incidence  (with  respect  to  the  target  normal).  The  total  neutron  yield  (last 
row  of  figures)  is  equal  to  the  neutron  yield  per  deuteron  multiplied  by  the  total  number  of 
deuterons  Ad+  =  tt  /?qLd^d  =  6.8  x  10^^  in  the  rear  layer  of  the  primary  target.  With  the  gold 
layer  thickness  increasing  the  neutron  yield  decreases.  This  can  be  expected  since  both  the 
mean  deuteron  energy  decreases  with  the  gold  layer  thickness  increasing  and  the  high-energy 
tail  of  lEDF  becomes  depleted. 

The  neutron  yield  versus  variation  of  the  thickness  of  the  deuterium  layer  of  the  primary 
target  is  more  subtle.  The  thinner  deuterium  layer  favors  the  production  of  higher  energy 
deuterons  (figure  8,  first  row)  with  narrower  angular  distribution  (figure  8,  second  row).  With 
an  increase  in  the  deuterium  layer  thickness  the  neutron  yield  per  ion  declines  (figure  8,  third 
row);  however,  the  number  of  deuterons,  which  scales  as  Ad+  ~  Td,  increases  and  so  does 
the  total  neutron  yield.  The  trade-off  leads  to  a  broad  peak  with  an  optimum  thickness  of  the 
deuterium  layer  between  50  and  100  nm.  According  to  figure  8  the  neutron  yield  is  insensitive 
with  respect  to  the  deuterium  layer  thickness. 

The  next  sequence  of  figures  illustrates  the  neutron  yield  versus  the  input  laser  parameters. 
The  deuteron  energy  and  angular  distribution  function  for  four  peak  laser  intensities  varying 
from  Iq  =  10^^  W  m~^  to  /q  =  10^^  W  m~^  are  plotted  in  figure  9.  We  first  make  a  comparison 
of  the  simplistic  model  with  the  Monte  Carlo  simulations  in  order  to  explore  the  role  of  the 
lEDF  on  the  neutron  production.  At  peak  laser  intensity  of  Iq  ~  10^^  W  m~^  the  neutron 
yield  from  the  Monte  Carlo  code  is  much  higher  compared  with  that  of  the  simplistic  model, 
while  at  peak  laser  intensity  of  /q  ~  10^^  W  m~^  the  exact  opposite  is  observed:  both  models 
yield  essentially  the  same  result.  This  behavior  can  be  explained  by  comparing  the  lEDFs 
for  both  intensities.  For  the  higher  peak  laser  intensity  Iq  ~  10^^  Wm“^  the  ions  form  a 
beam  with  a  relatively  broad  energy  spectrum  (figure  9(d)),  which  is  in  contrast  to  the  two- 
temperature  Maxwellian  lEDF  at  the  lower  peak  intensity  (figures  9(a)-(c)).  The  comparison 
clearly  shows  the  role  of  the  lEDF  on  the  neutron  yield.  The  most  remarkable  feature  of 
the  angular  distribution  of  the  neutron  yield  is  the  location  of  the  most  favorable  direction  for 
neutron  production.  At  low  laser  intensity  it  has  an  angular  dependence  close  to  a  sine  function, 
peaking  around  60-120°,  indicative  of  uniform  distribution  in  space  (figures  9(m)-(n)),  but 
with  an  increase  in  the  peak  laser  intensity  the  peak  shifts  to  smaller  angles  (figure  9(o))  and 
for  the  highest  laser  intensity  it  peaks  in  the  direction  of  laser  propagation  (figure  9(p)).  The 
increase  in  the  laser  pulse  duration  has  an  impact  similar  to  that  of  increasing  the  peak  laser 
intensity.  Both  the  neutron  yield  increase  and  the  angular  distribution  become  more  forward 
directed  with  the  laser  pulse  duration  increasing  (figure  10). 

A  summary  is  given  in  figure  11,  where  the  average  incident  deuteron  energy 
average  penetration  depth,  the  total  neutron  yield  per  ion,  the  number  of  deuteron  ions  and  the 
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Figure  8.  Deuteron  energy  distribution  function  (a)-(d),  angular  distribution  function  (e)-(h), 
neutron  yield  per  ion  versus  scattering  angle  (i)-(l)  and  total  neutron  yield  versus  scattering 
angle  (m)-(p)  for  different  thicknesses  Ld  of  the  D  layer.  The  laser  and  target  parameters  are 
Iq  =  10^^Wm“^,  Xq  =  1 /xm,  Rq  =  3 /xm,  tq  =  80fs,  ^laser  =  1.13 J,  Lau  =  0.25 /xm, 
Uq  =  lOOric. 

total  neutron  yield  are  plotted  versus  thickness  of  the  gold  and  deuterium  layers,  peak  laser 
intensity  and  laser  pulse  duration.  For  the  intermediate  peak  laser  intensity  of /q  =  lO^'^Wm"^ 
the  neutron  yield  is  10^-10^  neutrons,  depending  on  the  laser  pulse  duration.  Our  simulations 
are  in  close  agreement  with  the  simulations  in  [17],  which  report  neutron  yield  of  1.4  X  10^  at 
comparable  conditions  (peak  laser  intensity,  laser  wavelength,  spot  size  and  plasma  density). 
The  most  noteworthy  characteristics  regarding  the  total  neutron  yield  is  the  sharp  increase  in 
magnitude  with  the  peak  laser  intensity  (figure  11(^)).  This  result  is  in  agreement  with  the 
simulations  performed  by  Toupin  et  al  [17],  who  observed  an  increase  in  the  neutron  yield  of 
four  orders  of  magnitude  when  the  peak  laser  intensity  increases  from  Iq  =  10^^  W  m~^  to 
/o  =  2  X  10^^  W  m~^  (cf  figure  1 1(^)).  Measurements  of  the  neutron  yield  in  the  same  interval 
of  peak  laser  intensities  show  a  similar  trend  [6]. 

The  parametric  study  demonstrates  a  direct  correlation  between  the  neutron  yield  and 
directionality  and  the  ion  beam  characteristics.  The  latter  are  determined  by  various  input 
parameters  (laser  and  target),  which  can  be  controlled,  but  they  are  also  subject  to  ‘external 
factors  ’ ,  not  accounted  for  in  the  present  model.  One  of  them  is  the  2D  Cartesian  geometry  in  the 
laser-target  deposition  model.  We  use  only  three  EM  field  components  (E^,  Ey  and  B^),  while 
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Figure  9.  Deuteron  energy  distribution  function  (a)-(d),  angular  distribution  function  (e)-(h), 
neutron  yield  per  ion  versus  scattering  angle  (/)-(/)  and  total  neutron  yield  versus  scattering  angle 
(m)-(p)  for  different  peak  laser  intensities  Iq.  The  laser  and  target  parameters  are  A.o  =  1  /xm, 
Rq  =  3  /xm,  To  =  80  fs,  Lau  =  0.25  /xm,  Ld  =  0.05  /xm,  =  lOOric. 


a  truly  3D  electromagnetic  model  follows  all  six  components.  Thus  three-dimensional  effects, 
such  as  current  filamentation  and  quasistatic  magnetic  field  generation  [30, 31],  and  angular 
distributions  of  the  accelerated  electrons  [32]  are  not  fully  accounted  for.  For  example,  a 
strong  self-generated  longitudinal  magnetic  field  has  been  observed  in  simulations  [31],  which 
would  confine  the  electrons  and  ions  and  reduce  the  divergence  of  the  deuteron  beam.  Another 
simplification  is  the  neglect  of  collisions,  though  it  is  still  a  reasonable  approximation  at  the 
laser  intensities  and  particle  kinetic  energies  (~MeV)  we  deal  with.  Since  the  deuterons 
already  form  a  fairly  collimated  beam,  none  of  the  above  simplifications  is  expected  to  affect 
the  final  results  on  neutron  generation. 

4.3.  Neutron  spectra 

There  are  a  variety  of  applications  for  the  spectrum  of  energetic  neutrons  ranging  from  nuclear 
cross  section  information,  neutron  imaging,  fusion  and  fusion  reactions  and  the  material 
reactions  that  produce  gamma  radiation.  Both  neutrons  and  neutron  activated  gamma  radiation 
are  finding  widespread  application  in  the  detection  of  explosive  and  nuclear  materials.  A  useful 
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Figure  10.  Deuteron  energy  distribution  function  (a)-(d),  angular  distribution  function  (e)-(h), 
neutron  yield  per  ion  versus  scattering  angle  (/)-(/)  and  total  neutron  yield  versus  scattering  angle 
(m )-(p)  for  different  laser  pulse  durations  To .  The  laser  and  target  parameters  are /q  =  10^^Wm“^, 
A.0  =  1  /xm,  Rq  =  3  /xm,  Lau  =  0-25  /xm,  Ld  =  0.05  /xm,  =  lOO^c. 


discussion  can  be  found  in  the  paper  by  Ledingham  et  al  [33].  Figure  12  shows  integrated 
neutron  spectra  for  peak  laser  intensities  Iq  =  10^^,  10^"^  and  10^^  Wni“^.  Neutrons  are 
emitted  with  energies 

E„  =  {y/2  +  19.6/£i  (MeV)  +  cos^  6»„  +  cos  0„)^£’i/8,  (7) 

which  depend  on  the  incident  deuteron  energy  [17].  At  low  the  neutron  energy  is  not 
sensitive  with  respect  to  the  incident  deuteron  energy  and  peaks  near  2.45  MeV,  while  for  high 
the  neutron  is  emitted  with  energy  comparable  to  that  of  the  incident  deuteron.  For  the 
lowest  laser  intensity  the  neutron  spectrum  is  centered  around  2.45  MeV.  This  is  because  the 
deuteron  incident  energies  are  well  below  1  MeV.  For  the  intermediate  intensity  the  deuteron 
energies  extend  to  "^10  MeV,  showing  a  double  peak  at  around  2  and  3.5  MeV.  For  the  highest 
peak  laser  intensity  of  /q  =  10^^  W  m~^  the  ions  have  incident  energies  between  2  and  50  MeV, 
resulting  in  a  neutron  spectrum  peaked  at  around  2.5  MeV  with  a  long  tail  extending  to  tens  of 
megaelectronvolts.  These  high-energy  neutrons  can  be  more  efficient  for  neutron  activation. 
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Figure  11.  Average  incident  deuteron  energy  (a)-(d);  average  penetration  depth  of  deuterons  in  the 
secondary  target  (e)-(h)  neutron  yield  per  ion  (/)-(/),  total  number  of  deuterons  Nj)+  =  Tr/^gLonD 
(m)-(p)  and  total  neutron  yield  (q)-(t)  for  different  Au  layer  thickness  (first  column)  D  layer 
thickness  (second  column),  peak  laser  intensity  (third  column)  and  laser  pulse  duration  (right 
column).  Symbols — simulations  by  Toupin  et  al  [17]. 
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Figure  12.  Neutron  spectrum  for  peak  laser  intensities  Iq  =  10^^  (a),  10^^  (b)  and  10^^  W  m~^  (c). 
Laser  and  target  parameters:  A.o  =  1  Ro  =  ^  /^ni,  tq  =  80  fs,  Lau  =  0.25  /xm,  Ld  =  0.05  /xm, 
Uq  =  lOOtic. 
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4.4.  Neutron  activation 

Our  final  goal  is  to  estimate  the  number  of  neutrons  hitting  the  test  target.  Due  to  the  wide 
variety  of  targets  (particularly  size)  and  distance  to  target,  we  will  limit  the  study  to  one 
typical  case.  We  selected  a  target  the  size  of  a  basketball  (radius  /^target  =  0.12  m)  located 
^target  =  10  m  away  from  the  neutron  source  (figure  1  (/?)).  For  typical  laser  and  primary 
target  parameters  (/q  =  lO^'^Wm"^,  ko  =  1 /xm,  Rq  =  3 /xm,  tq  =  80  fs,  F^iaser  =  1.13  J, 
Lau  =  0.25  /xm  and  Ld  =  0.05  /xm)  the  number  of  neutrons  on  the  test  target  is  approximately 

^target  ~  ^target  Fq,  (8) 

where  ^target  =  ^target /^target  =  0.012  rad  (~0.7°)  is  the  half-angle  from  which  the  target  is 
visible  from  the  neutron  source.  Inserting  the  value  of  the  neutron  yield  in  the  first  angular  bin 
Yq  =  4.4  X  10^  neutrons  rad“^  (figure  l{o)),  we  get  Ftarget  ~  5  x  10^  neutrons. 

5.  Conclusion 

Numerical  simulations  of  the  angular  distribution  of  neutrons  from  high-intensity  laser-thin- 
foil  interaction  using  a  3D  Monte  Carlo  code  are  presented.  The  study  covers  variations 
of  target  thickness,  peak  laser  intensity  and  laser  pulse  duration  of  the  primary  target.  The 
neutron  yield  per  ion  and  the  total  neutron  yield  are  calculated  as  a  function  of  these  parameters. 
The  results  are  compared  with  a  simple  model  of  neutron  production  from  a  mono-energetic 
deuteron  beam  normally  incident  on  a  target.  Our  findings  can  be  summarized  as  follows: 

(i)  With  increasing  peak  laser  intensity  or  laser  pulse  duration  the  magnitude  of  the  neutron 
yield  increases  nonlinearly.  The  neutron  yield  decreases  with  the  thickness  of  the  gold 
layer  increasing  and  it  is  insensitive  to  the  thickness  of  the  deuterium  layer  of  the  primary 
target. 

(ii)  The  angular  distribution  of  neutrons  shows  a  weak  dependence  on  the  target  thickness,  but 
the  forward  scattering  improves  significantly  with  increasing  either  the  peak  laser  intensity 
or  the  laser  pulse  duration.  The  optimum  conditions  for  forward  neutron  scattering  are  a 
thin  primary  target  (gold  layer  thickness  Lau  =  0. 1-0.5  /xm,  deuterium  layer  thickness 
Ld  =  0.01-0.1  /xm)  and  laser  intensity  /q  >  10^"^  W m“^. 

(ii)  The  typical  neutron  yield  is  10“^-10“^  neutrons/ion  and  10^-10^  neutrons  laser 
energy. 

(iii)  The  simulation  results  show  that  for  a  typical  arrangement  of  laser  and  target  parameters 
some  ~10^  neutrons  can  be  put  on  a  target  ^^10  m  away. 
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